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Abstract 

This  research  program  investigates  matrix  crack  initiation  and  subsequent 
propagation  in  fiber-reinforced  ceramic  materials  for  use  in  high-temperature  structural 
applications. 

Though  it  presents  a  formidable  manufacturing  challenge,  the  inclusion  of 
ceramic  fibers  promises  to  increase  fracture  toughness  and  improve  failure  modes 
through  crack  deflection,  fi’acture  bridging,  and  frictional  interface  slip.  Experimental 
observations  show  that  ceramic  composites  initially  fail  at  several  points  in  the  matrix 
and  along  the  interfaces.  These  small  cracks  and  inherent  processing  flaws  propagate  and 
coalesce,  forming  large  cracks  that  lead  to  component  failure.  Therefore,  an 
understanding  of  small  crack  growth  is  necessary  for  the  design  of  composite  systems 
which  delay  critical  crack  formation  and  which  fad  in  a  desirable  manner. 

The  Surface  Integral  and  Boundary  Element  Hybrid  (SIB  EH)  method,  supported 
by  experimental  observations,  has  been  developed  to  model  crack  growth  in  brittle  com¬ 
posite  systems.  The  surface  integral  method  models  fractures  as  a  piece-wise  continuous 
distribution  of  displacement  discontinuities.  When  combined  with  traditional  boundary 
element  methods,  the  technique  provides  an  efficient  tool  for  modeling  three-dimensional 
crack  growth. 

This  approach  has  been  used  to  model  matrix  crack  initiation  in  a  lithium 
alumino-silicate  (LAS)  glass-ceramic  that  has  been  reinforced  with  continuous  silicon 
carbide  fibers.  By  modeling  the  effects  of  crack  pinning  and  bridging,  interfacial 
debonding,  and  ftictional  interface  slip,  this  investigation  aims  to  determine  the  stresses 
required  for  matrix  crack  initiation  and  the  material  parameters  which  promote  ‘^cefuT 
failure  modes.  These  results  have  been  compared  to  existing  analytical  solutions  for 
small  crack  growth.  Results  of  this  investigation  are  expected  to  be  useful  in  developing 
guidelines  for  the  manufacture  and  design  of  ceramic  materials  for  high-temperature 
structural  applications. 
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Chapter  1 
Introduction 


This  investigation  provides  a  fully  three-dimensional  analysis  of  small  matrix 
crack  growth  in  brittle  materials  reinforced  with  continuous  brittle  fibers  for  use  in  high- 
temperature  structural  applications.  The  project  focuses  on  the  influence  of  the  fiber- 
matrix  interface  on  small  matrix  crack  growth  and  on  fiber  failure  and  is  intended  to 
supplement  existing  models  by  assessing  the  three-dimensional  and  bimaterial  effects. 
Although  the  approach  developed  is  general,  a  particular  composite  system  consisting  of 
lithium  alumino- silicate  (LAS)  reinforced  by  silicon  carbide  fibers  has  been  modeled  to 
facilitate  comparison  with  experimental  data  and  existing  models.  A  modified  analytical 
expression  is  suggested  for  crack  initiation  and  subsequent  propagation. 

1.1  Ceramics  as  Structural  Materials 

Ceramics  are  attractive  stmctural  materials  because  they  offer  high  specific 
strengths,  excellent  thermo-mechanical  properties,  chemical  and  environmental  stability, 
and  low  raw  material  costs  [1].  For  high-temperature  applications,  such  as  internal 
combustion  engines  for  automotive  and  aeronautical  propulsion,  the  use  of  ceramics 
offers  great  gains  in  efficiency  because  of  their  insulating  properties  and  low  thermal 
expansion  coefficients.  In  many  proposed  applications  (e.g.  aerospace  propulsion 
systems  and  alternate-cycle  nuclear  reactors),  ceramics  are  the  only  possible  material 
choice  due  to  extreme  operating  temperature  requirements.  In  spite  of  these  demands. 
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ceramic  materials  have  found  limited  use  in  critical  structural  applications  due  to  their 
inherent  brittle  failure  modes  and  notch  sensitivity. 

In  the  spirit  of  Griffith's  brittle  fracture  investigations,  research  efforts  aimed  at 
improving  the  reliability  of  structural  ceramics  have  shifted  in  focus  from  flaw  control 
(i.e.  minimizing  the  processing  flaw  size  and  density)  to  damage  tolerance  (i.e.  improving 
the  material's  resistance  to  existing  flaws).  The  recognition  of  resistance  curve  behavior 
in  zirconia  and  the  development  of  high-temperature  reinforcing  fibers  sparked  this  shift 
and  has  renewed  interest  in  structural  ceramics  designed  for  toughened  fracture  behavior. 
Subsequent  experimental  investigations  and  theoretical  models  have  identified  the 
primary  toughening  mechanisms  and  quantified  many  of  their  effects  [2]. 

These  dominant  toughening  mechanisms  are  depicted  in  Figure  1.1  and  include 
transformation  toughening,  microcracking,  and  reinforcement  by  ductile  or  brittle 
inclusions.  The  first  two  are  process-zone  mechanisms  which  contribute  stress-induced, 
volumetric  dilation  or  softening  (respectively)  in  the  region  of  the  crack  tip  and  act  to 
shield  the  flaw  from  the  imposed  global  stress  state.  These  mechanisms  contribute 
modest  gains  in  overall  material  toughness  [2,3]. 

Additional  gains  in  toughness  can  be  achieved  through  bridging  mechanisms, 
which  transfer  a  portion  of  the  applied  load  to  intact  inclusions  spanning  the  fracture 
opening.  This  load  'shedding'  acts  to  reduce  the  stress  concentrations  at  the  crack  tip  and 
is  accomplished  through  shear  stresses  acting  across  bonded  or  frictionally  constrained 
fiber/matrix  interfaces.  Ductile  bridging  of  the  fracture  by  bonded  metallic  particles  can 
enhance  toughness  significantly  but  is  limited  to  lower  temperature  ranges  by  the  melting 
point  and  chemical  reactivity  of  the  inclusion  [3].  Whisker  or  particle  reinforcement  by 
brittle  inclusions  also  results  in  improvements  in  material  toughness  (roughly  200-300%) 
and  offers  higher  operating  temperature  ranges  [2]. 
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Figure  1.1  Dominant  Toughening  Mechanisms  for  Ceramic  Materials;  Dominant 
toughening  mechanisms  for  &ittle  ceramic  materials  include  (a)  transformation  toughen¬ 
ing,  (b)  microcracking,  (c)  ductile  reinforcement,  (d)  brittle  whisker  reinforcement,  and 
(e)  brittle  fiber  reinforcement. 


The  most  significant  increases  in  toughness  can  approach  two  orders  of  magnitude 
and  are  realized  through  continuous  fiber  reinforcement.  Although  it  presents  formidable 
manufacturing  challenges,  the  incorporation  of  brittle  fibers  also  offers  material  designers 
opportunities  for  damage  tolerant  behavior  and  more  'graceful'  failure  modes.  In  addition 
to  the  stress  relief  described  above,  local  nonlinear  effects  such  as  debonding  and  fiic- 
tional  slip  along  fiber/matrix  interfaces  can  contribute  to  the  overall  material  toughness 
[1,2].  The  strong  influence  of  interfacial  properties  on  these  toughening  mechanisms  has 
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been  demonstrated  and  suggests  tremendous  opportunities  for  designing  ceramic 
materials.  However,  tailoring  of  the  interface  for  toughened  behavior  remains  the  single 
greatest  challenge  for  successful  development  of  ceramic  composites. 

1.2  Failure  of  Fiber-Reinforced  Ceramics 

1.2.1  Damage  Development 

While  little  toughening  occurs  during  the  initial  stages  of  damage  development, 
this  process  has  a  significant  impact  on  the  final  failure  mode  and  on  the  toughening 
mechanisms  available.  Experimental  observations  show  that  'large'  matrix  cracks  (i.e. 
fractures  that  span  many  fibers)  initiate  from  manufacturing  flaws  in  the  matrix  that  are 
typically  on  the  order  of  the  fiber  spacing  [4,5].  These  small  cracks  propagate  under 
applied  loads  and  coalesce,  forming  large  fractures  which  eventually  lead  to  component 
failure  [6,7]. 

Theoretical  investigations  of  fracture  growth  in  fiber-reinforced  ceramics  in  which 
the  fiber  and  matrix  are  firmly  bonded  have  shown  that  the  stresses  imposed  on  the  fiber 
by  the  crack  tip  are  magnified  significantly  during  small  crack  growth  [8-10].  These 
local  stresses  pose  a  risk  to  the  integrity  of  the  fibers  and  therefore  to  the  toughened 
failure  modes  which  rely  on  them.  To  isolate  the  fibers  from  this  stress  concentration  and 
avoid  early  fiber  failure,  the  interface  must  be  sufficiently  weak.  Unfortunately,  this 
isolation  also  reduces  the  effectiveness  of  toughening  and  crack  pinning  mechanisms  [7]. 

For  these  reasons,  an  understanding  of  the  behavior  of  small  crack  growth  in 
ceramic  materials  reinforced  with  fnctionally-constrained  fibers  is  required  to  balance 
these  competing  effects  and  to  properly  tailor  the  interface  for  toughened  behavior.  The 
results  of  this  investigation  are  expected  to  provide  stractural  engineers  with  the 
necessary  tools  to  evaluate  the  integrity  of  ceramic  composite  components  and  to  provide 
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material  scientists  with  an  understanding  of  how  interfacial  effects  affect  the  overall 
composite  behavior. 

1.2.2  Failure  Modes 

The  ultimate  failure  of  composite  material  systems  is  complicated  by  fiber 
reinforcement,  but  has  been  categorized  by  Luh  and  Evans  according  to  the  macroscopic 
fracture  behavior.  These  failure  mode  classes  (Depicted  in  Figure  1.2)  can  be 
distinguished  by  the  progression  of  constituent  damage  and  are  dependent  on  the  relative 
structural  properties  -  predominantly  the  fiber  strength  and  interfacial  shear  strength  [7]. 

When  the  failure  strain  of  the  fibers  is  less  than  that  of  the  matrix  or  when  defects 
are  induced  in  the  fibers  during  processing,  the  failure  mode  is  dominated  by  frictional 
pullout  of  the  remaining  intact  fiber  segments  under  interfacial  shear  stress.  In  this  case, 
the  composite  fails  in  a  manner  similar  to  ceramics  reinforced  with  high  aspect-ratio 
whiskers  and  the  toughness  increase  AK  is  governed  by  the  interfacial  shear  stress,  X,  and 
the  average  pullout  length,  1.  Even  when  the  matrix  fails  preferentially,  fiber  failure  can 
occur  if  the  load  transferred  to  the  fibers  as  the  matrix  fractures  exceeds  the  fiber  strength. 
In  this  case,  progression  of  the  dominant  matrix  fracture  is  followed  closely  by  fiber 
failure  as  the  strong  interfacial  bonding  transfers  excessive  loads  from  the  failed  matrix  to 
bridging  fibers  in  the  crack-tip  wake.  While  both  of  these  modes  can  exhibit  increased 
toughness,  the  composite  remains  notch  sensitive  and  fails  in  catastrophic  manners. 

Only  when  the  failure  strain  of  the  matrix  is  lower  than  that  of  the  fibers  and  the 
fibers  are  sufficiently  strong  to  support  the  stresses  transferred  from  the  matrix  does  the 
ceramic  exhibit  notch  insensitivity  and  fail  in  a  'graceful'  manner.  In  this  case  -  depicted 
in  Figure  1.2(c)  -  damage  development  occurs  according  to  the  process  described  above 
in  Section  1.2.1,  leaving  a  matrix  riddled  with  cracks,  but  supported  by  intact  fibers. 
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Figure  1.2  Failure  Mode  Map  Failure  of  brittle  ceramics  reinforced  with  continuous 
brittle  fibers  is  complicated,  but  can  be  classified  according  to  the  macroscopic  fracture 
behavior.  These  failure  modes  are  dependent  on  the  constituent  microstructural  proper¬ 
ties  -  predominantly  the  ultimate  fiber  strength,  S,  and  the  interfacial  shear  stress,  x. 
(Adapted  from  Reference  [7]) 


More  importantly,  the  composite  may  continue  to  support  imposed  loads  in  spite  of 
extensive  damage  until  the  defects  can  be  detected  during  regularly  scheduled 
maintenance  operations. 

It  is  worth  noting  that  the  situation  depicted  in  Figure  1.2  is  further  complicated 
by  environmental  factors  (e.g.  oxidation,  irradiation),  material  changes  (e.g.  crystal 
growth  in  fibers),  and  stochastic  manufacturing  defects.  All  of  these  effects  can 
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significantly  influence  the  composite  behavior,  but  are  not  accounted  for  in  analytical 
fracture  models. 


1.2.3  Large  Crack  Models 


Experimental  evidence  and  theoretical  models  show  that  large’  matrix  cracks  (i.e. 
fractures  which  span  several  fibers)  bridged  by  intact  fibers  propagate  at  a  cracking  stress 
which  is  independent  of  crack  length.  Analytical  models  based  on  fracture  mechanics 
and  on  energy  considerations  have  been  developed  for  these  fractures  and  show  good 
agreement  with  experimental  data  [11-17].  Commonly  used  models  include  those 
derived  by  Aveston,  Cooper,  and  Kelly  (ACK),  by  Budiansky,  Hutchinson,  and  Evans 
(BHE),  and  by  Marshall,  Cox  and  Evans  (MCE)  [12,13,15].  These  theories  simulate  the 
dominant  toughening  mechanisms  which  occur  at  this  scale  -  fiber  bridging  and  frictional 
interface  slip  -  by  applying  uniform  distributed  closure  pressures  to  an  unbridged  fracture 
model. 

The  ACK  model  predicts  the  steady-state  matrix  cracking  stress,  as  a 

function  of  the  microstructural  composite  properties: 
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In  Equation  (1.1),  Ec  represents  the  composite  modulus  (EmVm+EfV^  calculated  using 
the  elastic  moduli,  E^  and  Ef,  and  volume  fraction,  Vm  and  Vf,  of  the  constituent 
materials.  The  fiber  radii,  R,  critical  interfacial  shear  stress,  x,  and  the  surface  energy  of 
the  matrix  material,  Ym,  also  influence  the  matrix  cracking  stress.  Similar  results  are 
obtained  using  the  BHE  and  MCE  models. 
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1.2.4  Small  Crack  Models 

Large  crack  models  have  been  extended  or  adapted  to  estimate  the  stress  required 
to  initiate  and  to  propagate  small  cracks  [15-17].  However,  as  the  fracture  size  decreases 
so  does  the  appropriateness  of  the  uniform  corrective  pressures  applied  in  these  models . 
In  addition,  small  fractures  may  be  subject  to  additional  toughening  mechanisms  not 
significant  for  larger  cracks  (e.g.  crack-tip  pinning  and  frictional  interface  slip  ahead  of 
the  fracture)  [2].  The  limits  of  applicability  for  these  'steady-state'  models  have  been 
estimated  to  be  several  fiber  spacings  or  greater.  Meda  and  Steif  have  demonstrated  the 
limited  applicability  for  these  fi’actiu’e  mechanics  models  when  the  radial  fracture 
dimensions  are  less  than  the  transition  flaw  size,  Ct  [18,19]: 
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In  Equation  (1.2)  Kjc  is  the  critical  mode  I  stress-intensity  factor  for  propagation  of 
bridged  matrix  flaws.  This  estimate  is  based  on  comparison  of  distributed  spring  models 
with  the  more  simplistic  long-crack  models  when  the  matrix  cracking  stress  is  within 
50%  of  the  ACK  cracking  stress.  The  transition  flaw  size  can  be  significantly  higher  than 
this  if  closer  agreement  is  required. 

The  acknowledged  disadvantage  of  these  models  is  that  they  do  not  accmately 
capture  the  toughening  mechanisms  that  occur  at  this  level,  including  the  stiffening 
effects  of  fibers  near  and  ahead  of  the  fracture  tip,  the  influence  of  interfacial  sliding  in 
this  region,  and  the  three-dimensional  nature  of  the  material  and  crack  propagation. 
Further,  these  models  make  no  estimate  of  the  load  transferred  to  the  fibers  or  of  the  local 
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Figure  1.3  Small  Matrix  Crack  Toughening  Mechanisms;  The  influence  of  crack  pin¬ 
ning  by  proximal  fibers  and  interfacial  slip  ahead  of  the  crack-tip  can  be  significant  for 
small  matrix  ftacture  propagation. 


stress  concentration  in  the  fibers  due  to  the  proximal  crack  tip  and  transferred  across  the 
frictional  interface. 

These  mechanisms  have  been  investigated  individually  in  detail,  beginning  with 
the  two-dimensional  work  by  Cook  and  Gordon  [20].  Subsequent  studies  have  captured 
new  aspects  of  the  problem,  including  line- tension  models  [21-23],  computational 
investigations  [8,10,24-27],  and  experimental  efforts  [22,28]. 

Several  analytical  models  have  been  developed  to  estimate  the  matrix  cracking 
stress  for  small  crack  growth  subject  to  these  toughening  effects  and  are  plotted  in 
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Figure  1.4  Toughened  Ceramic  Composite  Load-Deflection  Behavior;  Experimental 
observations  show  that  matrix  fracture  growth  begins  well  before  the  onset  of  visible 
nonlinearity  in  the  loda-deflection  curce  [4]. 


normalized  form  in  Figure  1.5.  Marshall,  Cox,  and  Evans  derive  the  following  form 
based  on  a  distributed-spring  model  [15]: 
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Similar  results  have  been  obtained  by  McCartney  and  are  included  in  Figure  1.5  [16]. 

The  effects  of  these  models  have  been  simply  summarized  by  Spearing  and  Zok  in 


the  following  relation  [29]: 
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Meda  and  Steif  have  further  improved  upon  these  relations  for  intermediate  length 
fractures  and  bridged  the  gap  between  the  small  crack  and  the  steady-state  crack  models 
using  an  axisymmetric  model  with  volume-weighted  effective  fiber  properties  [18-19]. 
Their  theory  links  the  matrix  regions  on  opposing  edge  of  the  'fiber'  to  simulate  the 
connection  which  exists  in  three-dimensions.  Although  no  explicit  relation  is  given  for 
the  penny-shaped  fracture,  normalized  results  are  included  in  Figure  1.5. 


1 .4  Research  Scope 

This  investigation  aims  to  contribute  understanding  of  the  influence  of  three- 
dimensional  and  bimaterial  aspects  of  the  problem  on  small  matrix  crack  behavior  using 
computational  fracture  mechanics.  The  focus  of  this  investigation  is  accurate 
determination  of  the  matrix  cracking  stress  for  small  fractures  in  a  stiffened,  three- 
dimensional  environment  and  investigation  of  the  influence  of  the  interface  properties  on 
the  fiber  stresses.  Comparison  of  these  results  with  existing  solutions  will  provide  some 
understanding  of  the  three-dimensional  and  bimaterial  effects  which  have  been 
approximated  in  previous  investigations. 
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Figure  1.5  Matrix  Cracking  Stress  Estimates  for  Small  Fractures;  Results  of  numerical 
and  analytical  models  for  mamx  cracking  stress  in  brittle  composite  materials  are  plotted 
as  a  function  of  the  normalized  characteristic  radial  crack  dimension  (c/Cn).  The  stresses 
have  been  normalized  by  the  ACK  steady-state  cracking  value. 


Determination  of  the  composite  failure  properties  is  accomplished  by  first 
developing  an  efficient  computational  scheme  for  fracture  analysis  in  composite  media. 
This  technique  is  based  upon  the  surface  integral  method,  which  models  three- 
dimensional  fractures  in  infinite  media  and  was  initially  developed  for  hydraulic  fracture 
applications  [30].  The  surface  integral  method  is  combined  with  traditional  boundary 
element  methods  using  superposition  to  incorporate  the  effects  of  model  boundaries  and 
stiffening  fibers.  For  analysis  of  uniaxially  reinforced  composite  materials,  fracture  cells 
are  constructed  on  various  scales  to  accurately  model  both  crack  initiation  and  small 
crack  growth.  Results  generated  include  the  matrix-crack  stress  intensity  factors,  applied 
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load,  interfacial  sliding  work  and  areas,  and  the  elastic  input  energy.  The  results  of  this 
investigation  suggest  modifications  of  existing  crack  theories  for  small  crack  growth. 

1.4.1  Silicon  Carbide/  Lithium  Alumino-Silicate  System 

Although  the  approach  is  general,  a  particular  composite  system  consisting  of 
lithium  alumino-silicate  (LAS)  reinforced  with  silicon  carbide  fibers  will  be  modeled. 
This  particular  material  has  been  chosen  to  facilitate  comparison  of  results  with 
experimental  data  and  observations  [4,6,7].  Lithium  alumino-silicate  is  a  glass-ceramic 
consisting  of  small  crystal  particles  ranging  from  0.05  to  1  micrometers  in  size 
surrounded  by  residual  glass  phase.  The  material  is  particularly  well-suited  for 
application  to  ceramic  composites  because  of  its  low  porosity,  ease  of  formation,  and 
small  crystals.  LAS  fractures  primarily  by  grain  boundary  cleavage  or  by  separation  of 
grain  clusters  and  can  operate  at  temperatures  slightly  above  1000°  C  [3,7,31]. 

The  reinforcing  fibers  are  silicon  carbide  fibers  and  are  usually  formed  by 
deposition  of  a  reactive  ceramic  on  a  fine  tungsten  core.  The  result  is  a  fiber 
approximately  15-20  micrometers  in  diameter  with  low  second-order  crystallinity. 
Although  their  small  size  gives  them  flexibility  necessary  for  processing,  these  'tows'  can 
be  affected  by  exposure  to  high-temperatures  and  irradiation  because  of  the  manufactured 
microstructure  [3,31]. 

1.4.2  General  Assumptions 

Several  general  assumptions  have  been  made  for  this  investigation  involving 
aspects  of  the  model.  Both  the  matrix  and  fiber  materials  have  been  modeled  as  linear- 
elastic,  isotropic  materials.  Although  the  matrix  material  is  subject  to  limited  plastic 
effects,  including  creep  deformation  and  microcracking,  these  effects  are  assumed  to  be 
localized  so  that  linear  -elastic  fracture  mechanics  is  applicable.  Although  the  relative 
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signMcance  of  these  effects  will  increase  with  decreasing  model  scales,  this  investigation 
will  provide  a  basis  from  which  to  evaluate  the  effects.  The  small  particle  size  and 
amorphous  structiue  of  the  matrix  and  fiber  (respectively)  suggest  that  this  assumption  is 
valid,  though  it  may  not  be  true  for  ceramics  in  which  the  grain  size  approaches  the  fiber 
spacing  or  in  which  the  fiber  crystal  structure  is  well  defined. 

The  interface  between  the  fiber  and  matrix  has  been  modeled  -  for  simplicity  -  as 
having  no  thickness  and  subject  to  constant  shear  stress  sliding.  The  assumptions  are 
justified  in  that  the  interface  slip  captures  the  effects  of  the  interphase  and  that  any 
mechanical  effects  of  the  inteiphase  can  be  incorporated  into  the  fiber  model.  Although 
theoretical  analysis  has  demonstrated  the  important  influence  of  the  normal  stress  in 
interfacial  sliding  on  push-out  tests,  experimental  evidence  for  these  particular  materials 
suggests  that  the  constant  shear  stress  model  is  sufficient  since  the  interface  is  frequently 
under  tension  or  separated  slightly  [7,32-34].  In  addition,  the  effects  of  thermal 
expansions  can  be  significant,  but  have  not  been  included  in  this  initial  investigation  [7]. 
These  factors  may  be  included  in  the  model  during  future  investigations. 

1.4.3  Model  Development  and  Results  Presentation 

The  development  and  analysis  of  this  matrix  crack  model  are  presented  in  the 
following  chapters.  Chapters  2  and  3  outline  the  development  of  the  surface  integral  and 
boundary  element  hybrid  method  for  application  to  three-dimensional  fracture  analysis  in 
composite  media.  Chapter  2  outlines  the  fundamental  theory  used  to  derive  the  surface 
integral  method,  presents  this  fracture  model  in  its  current  formulation  and  demonstrates 
its  capabilities  for  modeling  three-dimensional  fractures  in  infinite  media.  A  general 
integration  scheme  and  error  estimator  utilized  in  this  investigation  are  presented  as  well. 

To  evaluate  the  influence  of  composite  fibers  on  small  matrix  crack  growth,  the 
surface  integral  method  presented  in  Chapter  2  is  combined  with  classical  boundary 
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element  methods  using  superposition.  This  chapter  develops  this  hybrid  technique  in  its 
current  formulation.  After  briefly  reviewing  the  fundamentals  of  the  boundary  element 
technique,  Chapter  3  emphasizes  the  modeling  features  relevant  to  composite  fracture 
mechanics,  including  boundary  conditions,  subregions,  and  interfacial  slip  zones. 

Further  details  of  the  underlying  theory  and  of  the  current  implementation  for 
both  computational  schemes  are  outlined  in  the  Appendices  A  and  B. 

The  application  of  the  computational  approach  to  matrix  crack  initiation  and 
subsequent  propagation  is  presented  in  Chapter  4.  Verification  of  the  routine  by 
comparison  to  small  crack  propagation  experiments  is  described,  followed  by  results  for  a 
fully  three-dimensional  analysis  of  small  crack  growth  in  silicon  carbide  fiber-reinforced 
lithium  aluiTuno-silicate  subject  to  remote  tensile  stresses.  Pertinent  results  are  compiled 
and  presented  along  with  an  assessment  of  the  SIB  EH  method  for  computational  fi'acture 
analysis. 

Suggestions  for  modification  of  existing  cracking  stress  models  are  included  in 
Chapter  5  along  with  a  discussion  of  the  implications  for  the  manufacture  and  design  of 
ceramic  composite  materials . 


Chapter  2 
The  Surface  Integral  Method 


Computational  fracture  mechanics  models  based  on  force  dipoles  or  displacement 
discontinuities  can  provide  accurate  and  efficient  crack  solutions  for  linear  elastic  mate¬ 
rials.  For  this  investigation,  one  such  technique  -  the  surface  integral  method  -  is  em¬ 
ployed  to  capture  the  three-dimensional  aspects  of  small  matrix  crack  growth.  This 
chapter  develops  this  technique  in  its  current  formulation  and  demonstrates  its  capabili¬ 
ties.  Additional  details  of  the  underlying  theory  and  its  implementation  for  matrix  crack 
analysis  are  described  in  Appendices  A  and  B  respectively. 

2.1  Surface  Integral  Fundamentals 

The  surface  integral  method  models  three-dimensional  fractures  in  linear  elastic 
materials  as  a  piece- wise  continuous  distribution  of  displacement  discontinuities.  This 
technique  derives  from  the  general  concept  that  local  material  phenomena  can  be 
efficiently  modeled  with  dipole  distributions  and  resembles  the  indirect  boundary  element 
analysis  in  formulation  [35,36].  Development  of  this  computational  scheme  has  been 
motivated  by  the  need  for  efficient  crack  growth  models  which  rely  on  accurate  fracture 
parameter  solutions  and  simplified  growth  logistics.  Although  originally  developed  for 
hydraulic  fracture  applications,  the  surface  integral  method  has  been  used  successfully  to 
model  arbitrary  two-  and  three-dimensional  crack  growth  in  engineering  materials,  shear 
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band  formation  in  granular  media,  and  interfacial  slip  in  composite  materials  [26,30,37- 
42]. 

The  governing  integral  equation  expresses  the  stress  state  in  the  material 
surrounding  the  crack  as  a  function  of  the  displacement  discontinuity  distribution: 

t(x)  =  nJ^y‘(x,05(0<lAj  (2.1) 

where  t(x)  represents  the  traction  components  at  some  point  x  with  normal  direction  n  in 
the  media  surrounding  the  fracture  surface.  Sc-  The  integrand  combines  the  crack-face 
displacement  distribution,  5,  and  a  fundamental  solution  'y®,  which  gives  the  stresses  due 
to  unit  opening  of  an  infinitesimal  tensile  or  shear  crack  [35]. 

The  influence  functions  7®,  on  which  the  method  is  based  are  derived  by 
differentiation  and  combination  of  elasticity  solutions  for  point  forces  acting  in  an  infinite 
homogeneous  medium  [43].  In  this  formulation,  an  infinitesimal  tensile  crack  is 
represented  by  a  combination  of  dipoles  as  shown  in  Figure  2.1(a).  The  dominant  dipole 
simulates  the  tensile  crack  opening,  whereas  the  additional  dipoles  are  included  to 
counteract  the  associated  Poisson  contraction.  A  corresponding  multipole  can  be 
constructed  for  the  infinitesimal  shear  crack  opening  and  is  depicted  in  Figure  2.1(b)  [35]. 

2.2  Discrete  Formulation 

For  many  practical  applications  an  analytical  representation  for  the  crack-face 
displacements  cannot  be  obtained.  Therefore,  the  exact  distribution  is  approximated  in  a 
piece- wise  manner  by  dividing  the  crack  surface  into  subregions  over  which  some  locally 
continuous  distribution  is  assumed.  As  in  classical  boundary  element  methods,  the 
estimated  local  distribution,  8®(C),  is  defined  by  the  crack-face  displacements  at  specific 
points  within  each  element,  6®,  and  shape  functions,  N®(Q  [40].  In  this  formulation  the 
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Figure  2.1  Tensile  and  shear  displacement  discontinuity  representation;  (a)  and  (b) 
depict  the  dipole  combinations  us^  to  represent  infinitesimal  tensile  and  shear  fractoe 
events  respectively.  In  (a),  the  dominant  tensile  dipole  is  supplemented  with  perpendicu¬ 
lar  dipoles  to  counteract  Poisson  contraction.  In  (b),  the  shear  dipoles  are  balanced  for 
moment  equilibrium.  Both  are  expressed  as  displacement  discontinuities  by  combination 
with  appropriate  material  parameters.  (Adapted  from  Reference  [35]) 


integral  relation  in  Equation  (2.1)  becomes  a  summation  of  integrals  taken  for  each 


elemental  region,  Se,  comprising  the  fracture  surface. 

»(x)  =  2.n£T‘(jt,C)8’(0dA( 

(2.2) 

where  5*(0  =  Xn“(05“ 

(2.3) 

To  determine  the  crack-face  displacement  distribution,  a  collocation  method  can 
be  employed  in  which  the  applied  boundary  conditions  are  enforced  at  a  distinct  number 
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of  points  (collocation  points)  on  the  crack  surface  [38.35].  This  results  in  a  linear  system 
of  equations  relating  the  crack-face  displacements  and  tractions: 

C„8,=t,  (2-4) 


where  C„  =  n£„r(x,.ON»>(C)dAj  (2.5) 

The  coefficient  matrix  terms  Cu  represent  the  traction  forces  t  at  collocation  point  I 
corresponding  to  unit  crack-face  displacements  5  at  collocation  point  J .  The  integration 
in  Equation  (2.5)  is  taken  over  the  elements  enclosing  point  J. 

The  linear  system  expressed  in  Equation  (2.4)  can  be  solved,  and  the  results  can 
be  combined  with  prescribed  shape  functions  to  obtain  the  approximate  crack  opening 
distribution.  Stresses  and  displacements  at  points  in  the  surrounding  media  can  then  be 
expressed  as  a  function  of  the  crack-face  displacement  distribution: 

t(x)  =  X  8'n|^„T‘(x,ON<’>(OdA5  (2.6) 


u(x)  =  X  8’nJ^„T'‘(x,0N'’’(0dA(  (2.7) 

where  V*  in  Equation  (2.7)  represents  the  fundamental  displacement  solution  giving  the 
displacements  at  a  point  x  in  the  elastic  medium  surrounding  an  infinitesimal  fracture 
event 

In  general,  the  integral  terms  in  Equations  (2.5-2.7)  can  be  handled  using  two- 
dimensional  Gaussian  quadrature.  However,  when  the  collocation  points  at  which  the 
tractions  and  prescribed  crack  openings  are  evaluated  coincide,  the  l/R^  singularity  of 
the  fundamental  stress  solution  makes  the  integral  intractable.  Even  for  cases  of  a 
proximal  sampling  point  (important  for  the  hybrid  concept  developed  in  Chapter  3), 
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purely  numerical  integration  with  reasonable  integration  orders  is  often  corrupted  with 
significant  computational  errors. 

For  planar  crack  elements  with  internal  collocation  points,  these  situations  can  be 
efficiently  handled  by  subtracting  an  integral  term  equivalent  to  a  rigid  body  motion  [35]: 

C„  =  (N<'>(0  -  N'’')dA  +  fiA 

=  "i.,  rCN'"®  -  N<'>)dA  -  rdA  +  N»>n fdA  (2.8) 

where  NqW  is  the  shape  function  value  at  the  singular  collocation  point  and  Sj  represents 
the  entire  fracture  plane. 

The  first  integral  in  Equation  (2.8)  is  now  defined  in  the  sense  of  a  Cauchy 
principal  value  and  can  be  computed  directly.  The  third  term  represents  a  rigid  body 
motion  of  the  entire  fi-acture  plane  and  contributes  a  finite  value  to  either  the  stress  state 
(no  contribution)  or  to  the  displacements  (No).  Evaluation  of  the  second  integral  term  is 
somewhat  more  complicated  but  can  be  accomplished  by  first  recognizing  that  the 
individual  terms  of  the  fundamental  solutions  are  products  of  radial  terms,  7fe(r,z),  and 
angular  terms,  Ye(6).  when  expressed  in  a  local  cylindrical  coordinate  frame.  In  this 
form,  the  radial  terms  can  be  integrated  analytically.  The  remaining  angular  integrals  are 
then  evaluated  with  one-dimensional  Gaussian  quadrature. 


Y(x,C)  =  SYe(0)YR(r,z) 

(2.9) 

ee 

J  Y(x,C)dA  =  J  Ye(6)  J  YR(r,z)rdid0 

(2.10) 

r(e) 


In  research  conducted  independently,  this  integration  approach  has  also  been  applied 
successfully  to  regularized  integrals  for  the  boundary  element  method  [44]. 
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Despite  these  complicated  integration  procedures  and  its  limitation  to  linear 
elasticity,  the  surface  integral  method  provides  several  advantages  over  conventional 
numerical  techniques.  Because  the  fundamental  equations  are  based  on  multipole 
solutions  (representing  infinitesimal  fracture  events),  the  surface  integral  technique 
accurately  captures  the  stress  singularities  near  the  crack  tip.  Crack-face  displacements 
and  stress  intensity  factors  can  be  determined  with  a  limited  number  of  low-order  crack 
elements  [40]. 

More  importantly,  only  the  fracture  surface  is  discretized,  which  reduces  the 
required  degrees  of  freedom  and  simplifies  crack  growth  logistics.  Extension  of  the 
fractiue  surface  to  simulate  crack  propagation  simply  involves  the  addition  of  elements 
and  periodic  surface  remeshing.  This  is  a  considerable  advantage  when  compared  to 
classical  finite  element  methods,  which  require  significant  mesh  refinement  and  frequent 
volumetric  remeshing  in  the  regions  surrounding  the  propagating  crack  tip. 

2.3  Current  Implementation 

The  technique  outlined  above  has  been  implemented  for  analysis  of  three- 
dimensional  fractures  in  linear  elastic  media.  Although  the  approach  is  valid  for  arbitrary 
crack  geometries  and  boundary  conditions,  the  model  has  been  simplified  for  this 
investigation  to  planar  matrix  flaws  subject  to  boundary  conditions  which  are  symmetric 
about  the  fracture  plane  (i.e.  tensile  crack  opening  only).  Extension  of  the  approach  to 
more  general  situations  is  straightforward  and  has  been  presented  [26,35]. 

To  approximate  the  crack  opening  distribution,  the  fracture  surface  is  subdivided 
into  elemental  regions  over  which  local  variation  forms  are  prescribed.  The  fracture 
analysis  code  currently  offers  a  variety  of  element  configurations  as  outlined  in  Figures 
2.2  and  2.3.  Each  fracture  surface  can  be  subdivided  into  three-  and  four-sided  elemental 
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Figure  2.2  Nine-Noded  Element  Geometry;  The  fracture  surface  is  subdivided  into 
three-  and  four-sided  elements  with  linear  and  parabolically  curved  boundaries.  Mapping 
from  the  element  reference  frame  (b)  is  based  on  the  bi-quadratic  Lagrange  functions. 


regions  bounded  by  straight  or  parabolically  curved  boundary  segments.  For  integration 
piuposes,  points  on  the  elemental  fracture  surface,  Xp,  can  be  related  to  the  local  element 
reference  frame,  4,  by  the  bi-quadratic  mapping  functions,  M“(^),  from  the  Lagrange 
family  [45]. 

=  (2.11) 

Of  course,  simpler  element  geometries  are  possible  but  are  included  as  specialized  forms 
of  this  basic  nine-noded  Lagrange  element. 

The  singularity  of  the  fundamental  stress  solution  precludes  an  isoparametric 
representation  for  the  crack  opening  distribution  except  in  very  rigorous  formulations 
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[10,36,46].  However,  experience  indicates  that  accurate  ftacture  parameter  and  crack- 
face  displacement  solutions  can  be  obtained  with  simpler  local  distributions  [26,35]. 
Three  basic  shape  functions  have  been  found  sufficient  and  are  summarized  in  Figure  2.3. 
Implemented  options  include  constant,  constant-linear,  and  special  crack-tip  distributions. 
The  approximate  local  crack  opening  is  then  given  as  a  function  of  the  crack-face 
displacements  at  internal  collocation  points,  6“,  and  the  prescribed  elemental  shape 
functions,  N“.  Although  the  use  of  internal  collocation  points  results  in  a  discontinuous 
crack  opening  at  the  element  boundaries,  moderate  mesh  refinement  significantly  limits 
the  extent  of  these  discontinuities.  In  fact,  this  incompatibility  serves  as  a  useful  error 
indicator  as  discussed  in  Section  2.5  below. 

Accurate  fracture  parameter  and  crack-face  displacement  solutions  rely  on  the  use 
of  special  crack-tip  elements  for  the  crack  periphery.  The  crack-opening  distribution 
within  these  elements  varies  as  a  function  of  the  distance  from  the  crack  tip  according  to 
the  first  two  terms  of  the  Williams  expansion,  and  [47].  From  this  assumed 

variation,  two  elemental  shape  functions,  N“(^),  can  be  derived  in  terms  of  displacement 
values  at  the  internal  collocation  points  [26], 

pfpr-pfpr 


pfpr-pfpr 


(2.12) 


where  p(0  represents  the  distance  from  the  crack  tip,  and  pi  and  p2  are  the 
corresponding  collocation  point  distances.  In  contrast  with  conventional  formulations, 
the  crack-tip  shape  functions  depend  on  the  actual  crack  tip  radius  and  only  indirectly  on 
the  local  element  coordinates. 
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(a)  Constant  Element  (b)  Linear  Element  (c)  Crack-Tip  Element 

S(^)  =  5^  Sa)  =  d^)/2+  S(^)  =  f(5\5^)pa2r  + 

Figure  2.3  Summary  of  Element  Shape  Functions;  Crack-face  displacement  distributions 
within  an  elemental  subregion  are  defined  as  a  function  of  the  crack  opening  at  internal 
collocation  points  and  the  prescribed  shape  functions.  Constant,  linear,  and  special  crack- 
tip  functions  can  efficiently  approximate  crack  opening  for  general  three-dimensional 
fracture  situations. 


More  importantly,  local  stress-intensity  factors  can  be  accurately  computed  from 
the  crack  opening  displacement  at  some  small  distance  fix)m  the  crack  tip,  po.  For  the 
tensile  crack  case  [35]: 


K  -  G 

'  (l-v)2V2po/7t 


(2.13) 
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2.4  Fracture  Model  Results 

To  demonstrate  the  capabilities  of  the  surface  integral  method,  crack-face 
displacement  and  stress-intensity  factor  solutions  are  presented  for  a  variety  of  three- 
dimensional  fracture  configurations.  Figures  2.4  and  2.5  show  crack  opening  profiles 
superimposed  on  typical  surface  discretizations  for  penny-shaped  and  elliptical  cracks 
(respectively)  subject  to  uniform  internal  pressure.  Table  2.1  shows  the  corresponding 
stress-intensity  factors  for  selected  points  along  the  elliptical  crack  periphery. 
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Figure  2.4  Penny-Shaped  Crack  Subject  to  Internal  Pressure;  Crack  surface  dis¬ 
cretization  and  crack-face  displacement  contours  are  shown  for  a  penny-shaped  crack 
model  subject  to  uniform  internal  pressure  [47]. 
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Figure  2.5  Elliptical  Crack  Subject  to  Internal  Pressure;  Crack  surface  discretization  and 
crack-face  displacement  contours  are  shown  for  an  elliptical  crack  model  subject  to 
uniform  internal  pressure.  Stress-intensity  factors  have  been  estimated  at  the  points 
shown  and  are  tabulated  below. 


Table  2.1  Stress-Intensity  Factors  for  Elliptical  Crack;  Calculated  stress-intensity  factors 
for  selected  points  along  the  crack  periphery  show  good  agreement  with  analytical 
solutions  [48]. 


Point: 

(Ki/Ko)calc. 

(Ki/Ko)theor. 

%  Error 

1 

0.799 

0.80 

0.1 

2 

0.831 

0.84 

1.1 

3 

0.868 

0.89 

2.5 

4 

0.913 

0.94 

2.9 

5 

0.927 

0.96 

3.4 
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The  robustness  of  both  the  method  and  the  integration  scheme  is  demonstrated  by 
the  remaining  two  fracture  models.  The  first  (Figure  2.6  and  2.7)  involves  two  collinear 
penny-shaped  cracks  and  demonstrates  the  method's  capabilities  for  simulating  multiple, 
interacting  fractures. 

For  the  second  model,  the  fundamental  solutions  for  a  crack  in  an  infinite 
homogeneous  medium  have  been  replaced  by  solutions  derived  for  fractures  along  a 
bimaterial  interface.  These  influence  functions  have  been  derived  from  elasticity 
solutions  for  a  point  force  in  one  of  two  dissimilar,  bonded,  semi-infinite  regions  [43]. 
Although  the  analytically  integrated  terms  were  more  complicated,  the  derivation  follows 
the  same  approach  used  for  the  homogeneous  case.  Solutions  for  a  wide  range  of 
material  combinations  show  good  agreement  with  analytical  solutions  (Figure  2.8). 

2.5  Numerical  Issues 

In  addition  to  the  integration  issues  addressed  above,  several  other  numerical 
artifacts  can  significantly  affect  solution  accuracy  and  convergence.  These  include  the 
collocation  point  distribution,  the  prescribed  local  variation  (e.g.  shape  function  order, 
surface  discretization),  and  the  applied  boundary  conditions.  In  general,  accuracy  will 
improve  with  increasing  collocation  point  density  and  uniformity,  increasing  shape 
function  order,  and  decreasing  variation  in  local  boundary  conditions.  For  efficient 
problem  solution,  crack  surface  discretization  must  be  tailored  to  capture  the  local  crack 
opening  distributions.  Meshing  heuristics  have  been  developed  regarding  relative 
element  sizes  and  types  by  comparing  a  range  of  representative  computational  solutions 
with  known  analytical  solutions. 

In  addition,  an  error  indicator  based  on  the  discontinuity  in  estimated  crack 
opening  at  the  element  boundaries  has  been  employed  here  to  supplement  these 
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Figure  2.6  Collinear  Penny-Shaped  Cracks  Subject  to  Internal  Pressure;  &ack  surface 
discretization  and  crack-face  displacement  contours  are  shown  for  two  collinear  penny¬ 
shaped  cracks  subject  to  uniform  internal  pressure.  Stress-intensity  factors  have  been 
estimated  at  the  points  shown  and  are  tabulated  below  [49]. 


Figure  2.7  Stress-Intensity  Factors  for  Collinear  Penny-Shaped  Cracks;  Calculated 
stress-intensity  factors  are  plotted  as  a  function  of  the  position  sdong  one  crack  periphery 
for  two  cases  (r/d  =  0.8  and  0.94). 
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predictive  model  definitions.  To  generate  displacement  distribution  plots,  the  crack 
opening  at  the  nodal  points  are  estimated  as  a  function  of  the  surrounding  elemental 
variations.  Since  these  estimates  can  differ  across  element  boundaries,  the  maximum 
difference  as  a  fraction  of  the  averaged  value  can  be  used  as  a  qualitative  error  indicator 
[24].  Small  relative  discontinuities  in  the  crack-face  displacements  indicate  a  sufficiently 
refined  model,  whereas  large  values  can  signal  difficulties.  Tests  indicate  that  for  stress- 
intensity  factor  estimates  with  less  than  5%  uncertainty,  the  local  relative  displacement 
discontinuity  should  be  no  larger  than  20%. 
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Figwe  2.8  Calculated  Crack  Opening  Displacements  for  Penny-Shaped  Interface  Flaws 
Subject  to  Internal  Pressure;  (normalized  with  respect  to  the  homogeneous  case)  Crack 
opening  solutions  for  a  penny-shaped  interface  flaw  subject  to  uniform  internal  pressure 
show  good  agreement  with  analytical  predictions  for  a  wide  range  of  material 
combinations  [50]. 


Chapter  3 

The  Surface  Integral  and  Boundary  Element 

Hybrid  (SIBEH)  Method 

To  evaluate  the  influence  of  composite  fibers  on  small  matrix  crack  growth,  the 
surface  integral  method  presented  in  Chapter  2  has  been  combined  with  classical 
boundary  element  methods  using  superposition.  This  chapter  develops  this  hybrid 
technique  in  its  current  formulation.  After  briefly  reviewing  the  fundamentals  of  the 
boundary  element  technique.  Chapter  3  emphasizes  the  modeling  features  relevant  to 
composite  fracture  mechanics,  including  boundary  conditions,  subregions,  and  interfacial 
slip  zones.  Further  details  of  the  underlying  theory  and  of  the  current  implementation  are 
outlined  in  the  appendices. 

3.1  SIBEH  Method  Fundamentals 

The  effectiveness  of  the  surface  integral  method  for  modeling  three-dimensional 
fractures  in  infinite  domains  has  been  demonstrated.  Previous  investigations  have 
combined  the  surface  integral  and  finite  element  methods  to  model  fractures  in  the 
presence  of  finite  component  boundaries,  symmetric  planes,  material  interfaces, 
contained  crack-tip  plasticity,  and  thermal  strains  [26,35,37-42].  These  analyses  have 
used  the  surface  integral  method  for  accurate  fracture  solutions  and  relied  on  the  coupled 
finite  element  models  to  account  for  component  and  material  effects.  Successful 
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applications  of  this  approach  include  hydraulic  fracture  for  oil  and  gas  recovery,  crack 
growth  in  engineering  materials,  and  thermo-elastic  fatigue. 

For  the  present  analysis,  the  surface  integral  method  has  been  combined  with  the 
boundary  element  method  in  a  similar  hybrid  formulation  (SIBEH).  The  matrix  crack 
(surface  integral)  model  has  been  superpQsed  with  boundary  element  models  of  the 
surrounding  matrix  and  proximal  fibers.  Although  it  results  in  fully-populated,  coupled 
coefficient  matrices,  this  formulation  avoids  the  complicated  volumetric  finite  element 
meshes  which  would  be  required  for  this  problem.  In  addition  to  the  fracture  surface, 
only  the  material  interfaces,  symmetric  planes,  and  loading  surfaces  are  discretized. 

3.1.1  Elastostatic  Boundary  Element  Model 

The  boundary  element  method  can  be  derived  as  a  'weak'  formulation  of  weighted 
residual  statement  for  linear  elastostatics  [51,52].  The  resulting  integral  equation 
(Somigliana's  identity)  expresses  the  displacements  at  a  point,  I,  in  the  modeled  region  as 
a  function  of  the  traction  and  displacement  distributions  along  the  region  bounds,  F. 

Ui  (x)  =  uVx,  Xp  )T(Xp  )dAr  -  p*  (x,  Xp  )U(Xp  )dAr  (3.1) 

where  u*  and  p*  represent  influence  functions  derived  from  elasticity  solutions  for  a 
point  force  in  an  infinite,  homogeneous  domain.  T(xp)  and  U(xp)  represent  the 
boundary  traction  and  displacement  distributions,  respectively. 

In  practice,  these  distributions  are  approximated  by  dividing  the  model  boundary 
into  distinct  elemental  subregions  over  which  some  low-order  distribution  behavior  is 
assumed.  In  this  way,  the  tractions  and  displacements  can  be  expressed  in  a  piece-wise 
continuous  fashion  in  terms  of  locally  based  distribution  (shape)  functions,  N(J),  and  of 
the  traction  and  displacement  values  at  specific  boundary  (collocation)  points,  Tj  and  Uj. 
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TrXp;  =  XN''Vxp}T, 


Wx,)  =  Xn'''(-x,>U, 

(3.2) 

Using  these  approximations.  Equation  (3.1)  becomes: 

Ui(x)  =  ^Tj^u7x,Xp)N^^Vxp)dAi- 

(3.3) 

Applying  Equation  (3.3)  to  each  boundary  collocation  point  results  in  a 

equations  which  can  be  used  to  solve  general  boundary  value  problems: 

linear  system  of 

Hi,U,=GbT, 

(3.4) 

Hu  = 

(3.5) 

(3.6) 

where  c(J)  contains  geometric  constants  dependent  on  the  local  boundary,  and  the 
integrations  are  taken  over  the  region  surrounding  the  collocation  point  J. 

Stresses  at  points  within  the  model  can  then  be  expressed  as  a  function  of  the 
applied  boundary  values,  the  associated  solutions,  and  the  derivative  kernel  functions  d* 
and  s*  [51].  The  traction  force  at  a  point  in  the  model  interior,  t(x),  for  a  specific  normal 
direction,  n,  is  given  by  the  Equation  (3.7): 

t(x)  =  XTjnJ^s*(x.Xp)N^^Vxp)dAr 

-XUinJrdVx.Xp)N''Vxp)dAr  (3.7) 

Extension  of  these  relations  to  include  body  forces,  thermal  strains,  and  contained 
plasticity  is  straightforward,  but  has  been  omitted  from  this  derivation  for  clarity. 


39 


Chapter  3:  SIB  EH  Method 


3.1.2  Hybrid  Formulation 

The  process  used  to  couple  the  surface  integral  and  boundary  element  models 
combines  the  fundamental  relations  of  each  technique  using  superposition  as  depicted  in 
Figure  3.1.  The  problem  of  a  finite,  fractured  body  under  applied  crack-face  and 
boundary  tractions  (Figure  3.1(a))  can  be  solved  directly  by  superposing  and  linking  the 
two  models  as  shown.  Although  the  model  discussed  involves  only  traction  boundary 
conditions,  extension  of  the  resulting  relations  to  mixed  boundary  value  problems  can  be 
accomplished  simply  by  partitioning  the  coefficient  matrices  and  rearranging  terms. 

The  surface  integral  method,  shown  in  Figures  3.1(b),  models  the  fracture  in  an 
infinite  homogeneous  domain,  whereas  the  boundary  element  model  in  Figure  3.1(c) 
handles  the  finite,  uncracked  component.  This  approach  uses  the  accurate  fracture 
modeling  capabilities  of  the  surface  integral  method  to  greatest  advantage  while  retaining 
the  generality  of  the  boundary  element  method.  The  surface  integral  equation  system  is 
are  constructed  as  before.  However,  corrective  tractions,  t®,  (evaluated  along  the  image 
of  the  fracture  in  the  boundary  element  model)  must  be  subtracted  from  the  applied 
tractions,  t,  to  ensure  satisfaction  of  the  overall  boundary  conditions. 

[CP}  =  {1}-{1=}  (3.8) 

These  corrective  tractions  can  be  expressed  in  terms  of  the  boundary  element 
displacements  and  tractions. 

{t‘'}  =  [D]{U'*}-[S]{T'*}  (3.9) 

where  [D]  and  [S]  represent  the  integrals  expressed  in  Equation  (3.7)  and  are  evaluated  at 
images  of  the  surface  integral  collocation  points  in  the  boundary  element  model. 
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Fracture  Model 

Boundary  tractions,  T 
Crack-face  tractions,  t 


T 

Boundary  Element  Model 
Boundary  tractions,  T 
Corrective  tractions,  T* 


+ 


Surface  Integral  Model 
Crack-face  tractions,  t 
Corrective  tractions,  t* 


Figure  3.1  Superposition  of  Surface  Integral  and  Boundary  Element  Models;  To  solve 
the  problem  of  a  finite,  cracked  body  subject  to  applied  crack-face  and  boundary  tractions 
(a),  the  surface  integral  method  fracture  model  (b)  and  the  bound^  element  method 
component  model  (c)  can  be  superposed  and  linked  to  ensure  satisfaction  of  applied 
boundary  conditions.  Mixed  boundary  value  problems  can  be  solved  as  well  simply  by 
•  partitioning  the  resulting  equation  system  and  rearranging  terms. 
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Similarly,  satisfaction  of  the  global  conditions  along  the  component  boundaries  is 
ensured  by  combining  the  applied  tractions  with  corrective  tractions  from  the  surface 
integral  model  for  points  interior  to  the  finite,  uncracked  model  in  Figure  3.1(c): 

{t’*}  =  {T}-{t'^}  (3.10) 

Using  the  integral  relations  in  Equation  (2.6),  these  corrective  tractions  can  be  expressed 
as  a  function  of  the  crack-face  displacements: 

{T"}  =  [!]{«}  (311) 

In  this  case,  the  terms  of  matrix  [J]  are  evaluated  along  the  images  of  the  boundary 
element  collocation  points  in  the  surface  integral  model. 

Equations  (3.10)  and  (3.11)  are  applied  to  the  boundary  element  integral  relations 
for  the  model  depicted  in  Figure  3.1(c): 

[H]{U“}-[G]{T'”} 

=  [G]{T}-[G]{T"} 

=  [G]{T}-[GIJ]{5}  (3.12) 

However,  displacements  along  the  boundary  of  the  original  problem  are  a  sum  of 
the  displacements  from  both  superposed  models,  a  distinction  which  is  critical  for  mixed 
boundary  value  problems. 

{U}  =  {U“}  +  {U'*}  (3.13) 

The  displacements  {Usi)  are  given  by  the  matrix  form  of  Equation  (2.7)  evaluated  at  the 
images  of  the  boundary  element  collocation  points  in  the  surface  integral  model: 


{U-}  =  [L]{5} 
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# 


Combining  Equations  (3.8)-(3.14)  gives  the  complete  hybrid  method  equation  system, 
relating  the  crack-face  and  boundary  displacements  to  the  applied  tractions. 


'C-DJ-l-SL  S 
GJ-HL  H 


4' 

[o  gJItJ 


(3.15) 


By  combining  the  fracture  modeling  capabilities  of  the  surface  integral  method  with  the 
versatility  of  the  boundary  element  method,  the  SIBEH  method  provides  an  efficient  and 
robust  tool  for  linear  elastic  fracture  mechanics.  The  technique  is  particularly  well  suited 
for  fracture  propagation  analysis  since  only  a  limited  number  of  terms  in  Equation  (3. 15) 
need  to  be  recomputed  as  the  crack  face  is  extended. 

3.1.3  Multiple  Region  Models 

To  incorporate  the  stiffening  effects  of  inclusions  (e.g.  fibers)  in  composite  media, 
the  hybrid  formulation  presented  above  can  be  extended  using  a  subregioning  approach 
common  to  boundary  element  analysis  [51].  Additional  regions,  either  cracked  or  intact, 
can  be  linked  to  the  main  model  by  enforcing  displacement  equality  and  traction 
continuity  across  the  interface.  For  boundary  values  at  corresponding  interfacial  nodes. 


U‘=U^ 

rpl  ^  _ <^2 


(3.16) 


where  the  superscripts  1  and  2  denote  the  two  bonded  subregions.  When  solved  directly, 
the  resulting  set  of  equations  takes  the  form: 
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where  the  partitions  of  each  of  the  subre^on  equation  systems  have  been  combined 
according  to  the  relations  in  Equation  (3.16).  Additional  subregions  can  be  combined  by 
relating  the  boundary  values  at  common  interfacial  nodes  in  a  similar  fashion.  Although 
this  system  of  equations  can  be  used  to  accurately  model  composite  media,  this  approach 
leads  to  a  large  proportion  of  zero  terms,  suggesting  the  existence  of  more  efficient 
solving  schemes. 

3.2  Current  Implementation 

This  hybrid  approach  has  been  implemented  for  analysis  of  small  matrix  crack 
growth  in  composite  media  with  the  capacity  to  model  the  fractured  matrix  and  up  to 
three  additional  particles  or  fibers.  To  establish  the  system  of  equations  expressed  above 
in  Equations  (3.15)  and  (3.17),  the  matrix  and  fiber  surfaces  are  divided  into  three-  and 
foiu’-sided  elemental  regions  bounded  by  linear  or  parabolicaUy-curved  boundaries. 

Points  on  the  elemental  boundary  surface,  Xp,  can  be  related  to  a  local  element  reference 
frame  using  the  bi-quadratic  mapping  functions  described  in  Chapter  2: 

=  (3.18) 

For  simplicity,  lower  order  elements  are  modeled  as  reduced  forms  of  these  nine-noded 
Lagrange  elements. 

The  singularities  of  the  primary  fundamental  solutions  u*  and  p*  are  weaker  than 
those  derived  for  the  surface  integral  method  by  one  order.  Because  of  this  difference. 
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singular  integrals  for  points  on  element  boundaries  can  be  defined  in  a  limiting  sense  and 
evaluated  with  specialized  integration  schemes.  Therefore,  the  bi-quadratic  mapping 
functions  are  also  used  to  represent  the  local  boundary  value  variations  in  Equations 
(3.2).  This  iso-parametric  formulation  permits  continuous  boundary  value  distributions 
and  reduces  the  collocation  point  density  required  for  accurate  solutions.  Use  of  the 
Lagrange  family  of  functions  leads  to  relatively  uniform  collocation  point  distributions 
and  therefore  more  accurate  solutions. 

3.2.1  Singular  Integration  Scheme 

Singular  integrals  occur  regularly  in  the  boundary  element  formulation  when  the 
source  point  and  the  element  over  which  the  fundamental  solutions  u*  and  p*  are 
evaluated  coincide.  Many  methods  for  evaluation  of  these  singular  terms  have  been 
proposed,  including  modified  quadrature  rules,  element  subdivision,  analytical 
representations,  and  coordinate  transformations  [54-59].  The  most  elegant  (and  most 
accurate)  integration  scheme  results  in  complete  regularization  of  the  integrand  and  can 
be  applied  to  general,  three-dimensional  elastostatic  elements  [44].  This  recently 
developed  approach  has  been  implemented  for  singular  and  near-singular  integrals  and  is 
summarized  here  (Additional  details  can  be  found  in  Appendix  B). 

Using  an  approach  similar  to  that  presented  for  the  surface  integral  method,  the 
singular  integrals  are  regularized  by  subtraction  of  a  first-order  Taylor  series  expansion  of 
the  integrand  about  the  collocation  point,  Q,  nearest  to  the  source  point,  P.  To  facilitate 
semi-analytic  integration  of  the  regularizing  terms,  this  integral  is  evaluated  in  a  plane 
which  is  tangent  to  the  element  at  Q. 

/V{:p,OL.,ff)dAj+ 
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(3.19) 


where  f(Co>C)  f(Co>C)  represent  either  of  the  fundamental  solutions  evaluated 
on  the  element  surface  and  on  the  tangent  plane  (respectively)  and, 


Lo(0  =  N(Co)  + 


dN 


aN 


(3.20) 


When  expressed  in  a  specific  local  coordinate  frame,  the  first  derivatives  in  Equation 
(3.20)  can  be  evaluated  as  simple  functions  of  the  mapping  function  derivatives  and  nodal 
point  coordinates.  Finally,  semi-analytical  integration  of  the  regularizing  terms  is 
accomplished  using  the  approach  outlined  above  for  fracture  solutions  in  which  each  term 
of  the  fundamental  solution  is  separated  into  radial  and  angular  components. 

Although  this  integration  scheme  permits  evaluation  of  stress  and  displacement 
fields  nearer  to  the  element  surface  than  previously  possible,  it  is  not  generally  the  most 
cost-effective  approach.  For  source  points  further  from  the  element  than  the  element 
dimension,  straightforward  two-dimensional  numerical  integration  gives  equivalent 
results  and  is  less  expensive  [44]. 

The  singular  integrals  associated  with  the  boundary  displacements  need  not  be 
evaluated  directly  when  considering  finite,  bounded  regions.  By  virtue  of  a  rigid  body 
displacement,  these  terms  can  be  equated  to  the  sum  of  all  other  elemental  integrals  for 
the  singular  source  point  [51]. 

Hn=-SHu 

i*l 


3.2.2  Boundary  Conditions 

In  general  three-dimensional  elastostatics,  three  independent  conditions  must  be 
specified  at  each  boundary  collocation  point.  These  are  typically  the  values  of  the 
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displacement  and/or  traction  components  in  the  global  reference  frame.  However,  several 
situations  occur  frequently  which  require  special  treatment,  including  Dirichlet  comers, 
symmetric  planes,  and  region  interfaces  (treated  in  Section  3.3  below). 

The  term  'Dirichlet  comer*  refers  to  a  point  along  the  model  boundary  at  which  the 
boundary  tractions  change  discontinuously.  In  three-dimensions,  any  edge  or  comer  can 
present  modeling  difficulties  since  more  than  two  unknowns  -  the  comer  displacements 
and  tractions  for  each  neighboring  face  -  must  be  handled  with  one  equation  system.  The 
analogous  situation  in  two-dimensions  is  depicted  in  Figure  3.2(a).  For  certain  combina¬ 
tions  of  applied  boundary  conditions,  the  existing  approach  is  sufficient.  However,  the 
problem  becomes  indeterminate  when  only  the  displacement  components  are  prescribed. 

Various  approaches  have  been  developed  for  this  situation,  including  nodal  point 
separation,  discontinuous  elements,  and  derivative  formulations.  In  this  implementation, 
distinct  nodal  points,  separated  by  a  small  distance,  are  entered  at  such  'comers'  to 
accommodate  this  difficulty.  As  shown  in  Figures  3.2(b)  and  3.2(c),  the  distance  between 
nodal  points  can  either  be  spanned  by  a  small  element  or  left  as  a  gap  in  the  component 
boundary.  Separation  of  the  collocation  points  reduces  the  problem  to  determinate  form 
by  increasing  the  number  of  available  equations.  However,  care  must  taken  to  separate 
the  points  sufficiently  to  create  distinguishable  equations  and  permit  accurate  integration, 
but  not  so  far  that  the  disraption  of  the  component  boundary  affects  the  model.  Typical 
values  used  for  the  matrix  and  fiber  models  range  from  1/4  to  1/2  local  element 
dimension. 

Planes  of  symmetry,  used  to  reduce  and  simplify  models,  can  present  similar 
difficulties.  Boundary  conditions  along  symmetric  planes  mandate  that  the  normal 
displacement  and  tangential  traction  components  be  zero.  When  these  planes  are  aligned 
with  the  reference  axes,  these  conditions  can  be  applied  as  prescribed  boundary 
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displacements  and  tractions.  However,  for  general  planes  of  symmetry,  the  boundary 
relations  become  more  complicated: 


U  n  =  0, 

T  •  r  =  0,  and  T  •  s  =  0 


where  r 


ro-(n-ro)n 

IK-(nToH 


and  s  =  rxn. 


(3.22) 

(3.23) 


In  Equation  (3.23)  Tq  is  some  initial  guess  for  the  tangent  vectors  r  and  s,  and  n  is  the 
outward  unit  normal.  These  additional  equations  are  combined  with  the  existing  system 
of  equations  to  relate  the  boundary  values  along  symmetric  planes  so  that  the  conditions 


//////////M  , 
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Figure  3.2  Treatment  of  Dirichlet  Comers;  Boundary  points  (a)  at  which  discontinuities 
in  the  applied  conditions  create  intractable  problems  (e.g.  comers  and  edges  subject  to 
applied  displacements)  can  be  handled  in  the  boundary  element  formulation  by  creating 
two  distinct  nodal  points,  joined  (b)  or  unjoined  (c),  at  the  point  of  discontinuity. 
Separation  of  these  nodes  should  be  sufficient  to  distinguish  the  resulting  equations  and 
to  permit  accurate  integration  but  should  not  dismpt  the  model.  (Adapted  from  [51]) 
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expressed  in  Equation  (3.22)  are  prescribed  in  the  solution. 

To  accurately  simulate  symmetric  boundaries,  corrective  tractions  applied  to  the 
boundary  element  model  must  reflect  the  variation  in  stress  due  to  the  proximal  fracture 
elements.  For  reasons  discussed  later  in  this  chapter,  the  local  accuracy  is  dependent  on 
the  refinement  of  the  boundary  mesh,  the  relative  surface  integral/boundary  element  size 
and  proximity  of  the  crack  elements  to  the  boundary.  To  reduce  computational  errors  for 
commonly  used  symmetry  conditions,  single-,  two-,  and  three-fold  symmetric  crack 
elements  are  used  for  near-boundary  fracture  elements.  These  symmetries  are 
incorporated  in  the  surface  integral  technique  by  reflecting  the  fracture  elements  about  the 
appropriate  planes  and  condensing  the  resulting  matrix  terms.  The  effect  in  the  hybrid 
method  is  to  translate  the  influence  of  the  fracture  on  the  neighboring  boundary  from 
shaiply-varying  traction  to  milder  displacement  distributions.  Because  the  evaluated 
terms  are  subject  to  computational  noise,  the  condensed  values  must  be  filtered.  Typical 
results  are  shown  in  Figures  3.5  and  3.6  for  the  case  of  a  centrally  located  penny-shape 
fracture  in  a  tensile  bar.  Even  for  low  degree-of-freedom  models,  use  of  symmetric 
fracture  elements  reduces  the  errors  in  stress-intensity  factors  (and  therefore  crack-face 
displacements)  for  a  wide  range  of  relative  element  sizes. 

3.3  Bimaterial  Interface  Models 

An  important  aspect  of  the  matrix  crack  model  is  the  effect  of  interfacial  slip  and 
debonding.  As  described  in  Chapter  1,  proper  tailoring  of  the  interfacial  properties  can  be 
critical  to  the  overall  composite  performance.  A  technique  borrowed  from  contact 
mechanics  analysis  has  been  employed  to  provide  a  flexible  means  of  incorporating  this 
phenomenon  in  the  SIBEH  method. 
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Figure  3.3  1/6  Symmetric  Model  of  a  Circular  Crack  in  a  Tensile  Rod;  The  surface 
integral  and  boundary  element  meshes  are  shown  for  the  1/6  symmetric  model  used  to 
evaluate  the  effectiveness  of  symmetric  fracture  elements. 


Figure  3.4  Symmetric  Model  Error  Estimates;  Effof  estimates  for  the  1/6  symmetric 
model  shown  in  Figure  3.3  are  plotted  as  a  function  of  relative  crack  radius.  In  this 
model,  the  relative  crack  radius  also  represents  the  relative  surface  integral/  boundary 
element  size,  important  for  accurate  fracture  parameter  solutions. 
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3.3.1  Interfacial  Slip  and  Separation 

The  interfacial  phenomenon  relevant  to  composite  damage  mechanics  can  be 
modeled  by  three  distinct  situations:  perfect  bonding,  interface  debonding,  and  frictional 
slip.  Bonding  occurs  during  manufacture  because  of  chemical  interaction  between  the 
fiber  and  matrix  or  due  to  thermal  strains  acting  across  the  interface  and  is  modeled  by 
enforcing  displacement  and  traction  continuity  as  described  above.  Interfacial  debonding 
or  sUp  are  simulated  by  applying  boundary  tractions  according  to  the  composite 
microstructure  and  the  local  conditions.  If  the  stresses  acting  normal  to  the  interface  are 
tensile  and  above  a  critical  value,  the  interface  debonds  and  the  applied  tractions  become 
zero.  This  critical  normal  stress.  On.  is  used  to  reduce  the  effects  of  small  numerical 
errors  in  interface  tractions  and  can  simulate  the  binding  influence  of  thermal  strains. 

The  bonding/  debonding  conditions  are  summarized  in  equation  (3.24): 


Bonding:  <7  <  <T„ 


Debonding:  <T  ><T„ . 


(3.24) 


Excessive  shear  in  conjunction  with  sub-critical  normal  stresses  leads  to 
interfacial  slip,  which  may  be  simulated  in  one  of  two  manners.  A  Coulombic  friction 
rule  is  commonly  used,  in  which  the  interfacial  shear  stress  is  proportional  to  the  local 
normal  stress  (according  to  friction  coefficient,  |x)  and  is  taken  in  the  direction  opposite 
impending  slip: 

T  for  T<Ii\o\  and  (T^(7„ 


x  =  - 


p\a\  for  T</t|a|  and 


(3.25) 


Interfacial  sliding  on  this  scale  for  fiber-matrices  is  small  so  that  it  may  be  appropriate  to 
use  the  static  or  break-away  friction  coefficient  rather  than  the  kinetic  coefficient.  These 
values  are  determined  from  fiber  push-out  or  pull-out  tests  [32-34]. 
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An  alternative  interface  model  is  constant  shear  stress,  in  which  the  critical  sliding 
stress,  XcRTT^  is  a  constant  model  parameter: 

'r=TcRrr  for  T<TcRrr  and  G<c^.  (3.26) 

As  in  the  Coulomb  friction  model,  the  sliding  shear  stress  acts  in  a  direction  opposite  slip. 
This  model  is  a  simpler  approximation  and  is  sufficient  for  many  ceramic  composite 
systems,  including  SiC/LAS  [34]. 

Although  the  iterations  required  for  convergence  may  differ,  the  approach  to 
solving  the  non-linear  interface  problem  is  the  same  for  both  interface  models.  The 
simplest  approach  is  to  completely  solve  the  non-linear  interfacial  slip  problem  for  each 
crack  advance.  For  each  fracture  configuration,  the  interfacial  nodes  are  sampled  and 
updated  according  to  the  criteria  expressed  in  Equation  (3.25)  or  (3.26)  until 
convergence.  For  crack  propagation  problems,  the  interfacial  slip  zone  is  recorded  and 
used  as  an  initial  condition  for  subsequent  fracture  geometries.  This  reduces  the  number 
of  iterations  required  for  small  fracture  steps,  and  improves  the  converged  solution. 

Even  with  special  matrix  solvers,  this  process  is  costly.  Improvements  in  the  non¬ 
linear  solution  strategy  can  lead  to  large  reductions  in  computational  effort.  Although  it 
has  not  been  implemented  here,  Larson  successfully  employed  an  iterative  scheme  which 
estimated  fracture  advance  and  interfacial  slip  alternately  [26].  Further  numerical  issues 
which  have  been  addressed  include  inteipenetration  of  opposing  faces  and  low-level 
oscillations  in  the  normal  stress. 

3.3.2  Experimental  Verification 

For  sufficiently  weak  interfaces,  interfacial  slip  actually  occurs  ahead  of  the 
fracture  tip  [2,26].  In  fact,  this  debonding  is  critical  to  reduction  of  fiber  stresses  for 
toughened  failure  behavior.  To  verify  the  accuracy  of  the  interfacial  model  employed,  an 
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interface  slip  zone  has  been  generated  and  observed  using  the  Interface  Slip  Experiment 
apparatus  developed  by  Larson  and  rebuilt  for  this  investigation  [26].  Experimental 
results  serve  as  a  qualitative  indicator  of  the  slip  zones  expected  on  curved  fiber-matrix 
interfaces  and  give  a  quantitative  assessment  of  the  computational  simulations. 

The  experimental  apparatus  depicted  in  Figure  3.5  consists  of  a  Sylastic  rubber 
block  constrained  between  two  'rigid*  PMMA  plates.  A  penny-shaped  flaw  approaching 
the  interface  is  simulated  with  an  inflatable  rubber  diaphragm  located  beneath  the  rubber 
block.  The  interface  conditions  can  be  varied  by  changing  the  confining  strains 
(measured  with  strain  gages)  and  the  interfacial  coating.  Although  the  rubber  material 
does  not  behave  like  a  brittle  ceramic,  its  use  allows  interfacial  deformations  to  be 
observed  directly  by  comparing  grids  marked  on  both  the  block  and  on  the  PMMA. 

Figure  3.6  shows  a  typical  experimental  slip  pattern  and  the  slip  zone  predicted 
using  the  1/4  symmetric  model  shown  in  Figure  3.7.  For  this  case,  the  interface  slip 
model  is  defined  by  a  constant  shear  stress  estimated  by  Larson  [26].  Both  the  predicted 
slip  zone  location  and  the  relative  displacements  show  good  agreement  with  experimental 
results.  Although  there  is  a  slight  difference  in  the  extent  of  the  slip  zone,  the  effect  of 
this  on  a  neighboring  fracture  is  expected  to  be  relatively  small. 


Figure  3.6  Experimental  and  Computational  Interface  Slip;  Predicted  slip  along  a  planar 
bimaterial  interface  show  good  agreement  in  general  shape,  location,  and  magnitude  with 
experimental  slip  zones  observed  using  the  Interface  Slip  Experiment  [26]. 
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Figure  3.7  1/4  Syrnmetric  Model  , of  Interface  Slip  Experiment;  Surface  integral  and 
boundary  element  subrfegion  meshes  are  shown  for  a  1/4  symmetric  model  of  the 
Interface  Slip  Experiment.  A  contour  plot  of  interfacial  slip  has  been  superimposed. 
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sensitive  to  the  collocation  distribution  (fracture  and  boundary  meshes)  and  -  to  a  lesser 
extent  -  the  integration  order.  Attempts  at  iterative  improvement  of  the  solution  resulted 
in  insignificant  changes,  indicating  that  the  solution  scheme  (LU  decomposition  with 
partial  pivoting)  performs  well,  even  for  large  problems  (>5000  degrees  of  freedom). 

This  initial  investigation  has  developed  and  made  use  of  several  meshing 
heuristics  to  improve  solution  accuracy.  These  guidelines  can  easily  be  incorporated  in 
mesh  generators,  although  they  have  not  yet  been  implemented  in  this  fashion.  Although 
some  flexibility  can  be  tolerated,  best  results  are  obtained  with  uniform  collocation  point 
density.  For  fracture  models,  this  means  relative  element  sizes  proportional  to  the 
number  of  collocation  points  per  element  as  discussed  in  Chapter  2.  Most  commonly, 
crack-tip  elements  should  be  roughly  twice  as  large  as  neighboring  constant  elements. 

Use  of  Lagrange  family  elements  and  uniform  mesh  densities  have  addressed  this  issue 
for  the  boundary  element  models.  Other  common  error  sources  include  severely  distorted 
elements  and  small  comer  angles  [45,51]. 

Additional  issues  arise  with  the  combination  of  the  siuface  integral  and  boundary 
element  methods.  The  fracture  surface  must  be  sufficiently  refined  to  accurately  model 
the  crack-opening  displacements.  The  boundary  element  mesh  should  reflect  local 
boundary  value  variations  due  to  the  presence  of  the  crack  and  the  applied  conditions. 
However,  the  integration  errors  for  both  techniques  increase  with  relative  source  point 
proximity  (See  Figure  3.4  above).  While  some  special  conditions  can  be  simulated 
directly  (e.g.  symmetric  fracture  elements),  these  conflicting  demands  must  generally  be 
balanced. 

3.4.2  Improved  Iterative  Scheme 

The  presence  of  a  significant  number  of  zero  terms  in  the  multi-region  equation 
system  described  by  Equation  (3.17)  suggests  opportunities  for  improved  storage  and 
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solution  schemes  [60,61].  For  this  investigation,  an  iterative  improvement  scheme  has 
been  developed  which  achieves  these  goals.  The  technique  applies  inteifacial  traction 
condition  only  to  the  larger  fracture-matrix  model  so  that  this  equation  system  can  be 
decomposed  and  stored.  Separate  systems  of  equations  are  setup  and  decomposed  for 
additional  subregions  (e.g.  fibers)  for  each  interface-slip  step.  Bonded  interface  points  in 
the  fiber  models  are  subject  to  the  corresponding  matrix  displacements.  The  routine 
iterates  between  the  models  in  this  fashion  until  convergence  to  guarantee  satisfaction  of 
the  interface  continuity  (See  Figure  3.10). 

For  the  cracked  matrix  model  at  solution  iteration  j: 


C-DJ-hSL  S‘ 
GJ-HL  H‘ 


(3.27) 


The  initial  interface  traction  conditions  are  taken  from  the  previous  fracture  analysis  in 
crack-growth  studies  to  reduce  number  of  iterations  required  for  convergence. 

The  corresponding  relations  for  the  attached  fiber  regions  can  be  expressed  in 
matrix  form  as: 


[H'  ={G^T^}+[-H""]{U""}, 


(3.28) 


The  iterative  scheme  is  effective  since  the  boundary  conditions  imposed  on  the  matrix 
interface  points  do  not  change  type  when  slip  occurs.  This  removes  a  decomposition  of 
the  main  equation  system  (N^  operation)  from  the  interface  slip  loop,  leaving  only  the 
substitution  required  for  each  solution  step  (N^  operations).  For  typical  matrix  crack 
model  sizes,  this  approach  results  in  a  reduction  of  effort  when  less  than  200  iterations  are 
required  for  convergence.  It  is  important  to  note  that  this  is  not  generally  true.  For 
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systems  with  stiffer  fibers,  the  interfacial  tractions  must  be  damped  to  guarantee 
convergence.  The  best  approach  depends  on  the  model  geometry  and  material  properties. 

3.4.3  Computational  Efficiency 

One  goal  of  this  research  program  has  been  to  develop  a  computational  fracture 
mechanics  scheme  which  can  be  operated  on  mid-range  computers.  While  the  SIBEH 
code  developed  for  the  matrix  crack  analyses  has  been  tested  and  run  on  CRAY 
supercomputers,  the  results  presented  in  this  document  were  generated  on  a  DEC-Station 
3100.  Large  degree-of-freedom  problems  have  been  solved  in  several  hours  with  the  use 
of  symmetric  condensations,  semi-analytic  integration  schemes,  and  the  iterative  multi¬ 
region  solution  scheme  described  above.  Since  most  of  the  computational  effort  involves 
integration,  storage  of  the  model  equation  systems  in  formatted,  direct-access  files 
significantly  reduces  the  time  required  for  repeated  analyses  of  the  same  model  subject  to 
different  boundary  or  interface  conditions  [62]. 
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The  application  of  the  computational  approach  developed  in  Chapters  2  and  3  to 
matrix  crack  initiation  and  subsequent  propagation  is  presented  in  this  chapter.  Results 
for  a  fully  three-dimensional  analysis  of  small  crack  growth  in  a  silicon  carbide  fiber- 
reinforced  glass-ceramic  subject  to  remote  tensile  stresses  are  outlined  with  particular 
attention  to  matrix  cracking  stresses  and  fiber  failure.  A  discussion  of  findings  and  their 
relation  to  other  matrix  crack  theories  is  included  as  well.  The  implications  for 
toughened  ceramic  materials  are  discussed  in  Chapter  5. 

4.1  Small  Matrix  Crack  Models 

The  composite  system  modeled  in  this  investigation  consists  of  a  glass-ceramic 
matrix  reinforced  by  40%  volume  silicon  carbide  fibers  aligned  with  the  loading  direction 
as  depicted  in  Figure  4.1.  The  constituent  material  properties  and  microstructural 
conditions  are  summarized  in  Table  4.1.  The  interface  is  modeled  with  no  thickness  and 
is  governed  by  constant  shear  stress  conditions  as  suggested  for  this  material  combination 
by  experimental  results.  Four  interfacial  strengths  will  be  evaluated,  ranging  from  lightly 
bonded  (2  MPa)  to  strongly  bonded  (40  MPa). 
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Figure  4.1  Matrix  Crack  Initiation  Configuration;  The  matrix  crack  model  simulates  an 
initially  circular  flaw  propagating  in  a  glass-ceramic  matrix  stiffened  by  hexagonally- 
packed  silicon  carbide  fibers.  This  particular  configuration  is  based  on  experimental 
observations  for  this  material  system  [4,5]. 


The  matrix  crack  is  assumed  to  begin  as  a  penny-shaped  flaw  located  between  the 
fibers  and  aligned  perpendicular  to  the  fiber  and  loading  direction.  This  configuration  is 
based  on  acoustical  and  optical  observations  of  crack  initiation  in  this  particular  material 
system  and  has  been  modeled  using  approximate  methods  to  determine  the  matrix 
cracking  stress  [4,5,19,27].  The  model  constracted  for  this  situation  describes  matrix 
crack  intiation  and  subsequent  growth  when  the  fracture  is  on  the  order  of  the  fiber 
spacing.  The  scope  has  been  chosen  to  bridge  the  gap  from  the  initiation  point  to  the 
axisymmetric  models  of  'intermediate'  size  (i.e.  several  fiber  spacings  in  size,  but  less 
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Table  4.1  SiC/LAS  Composite  Material  Properties;  [7,15] 


SiC  fiber: 

Ef 

200  GPa 

Vf 

0.2 

S 

2  GPa  at  room  temp. 

1  GPa  at  1000  C 

R 

8  |i,m 

LAS  matrix: 

Em 

85  GPa 

Vm 

0.3 

K,c«" 

2  MPaVm 

Ym 

21  J/m2 

Composite: 

Vf 

0.40 

X 

2  -  40  MPa 

than  critical  transition  flaw  size)  [18].  The  model  is  further  reduced  using  symmetric 
boundary  conditions  and  condensation  to  include  the  regions  depicted  in  Figure  4.2. 

4.1.1  Matrix  Crack  Initiation 

The  primary  focus  of  this  investigation  is  the  behavior  of  matrix  cracks  initiating 
between  fibers  and  propagating  past  the  first  rows  of  fibers.  The  1/6  symmetric  model 
incorporates  the  stiffening  effects  of  the  first  two  fibers  and  has  been  tested  for  three 
distinct  crack  lengths  subject  to  four  interface  conditions.  The  mesh  and  model  cell  are 
shown  in  Figure  4.3.  The  model  consists  of  roughly  2500  degrees  of  freedom  (including 
coincident  interface  points)  before  condensation  for  vertical  symmetries  and  requires 
slightly  more  than  2  hours  for  complete  interface  slip  analysis. 
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Figure  4.2  Small  Matrix  Crack  Configuration  and  Model  Cell;  Model  cell  Iwundaries 
are  constructed  using  symmetry  as  shown  for  investigation  of  matrix  crack  initiation  and 
subsequent  small  matrix  crack  growth. 


To  simulate  remote  applied  strains,  the  upper  and  lower  surfaces  are  constrained. 
More  accurate  results  are  obtained  for  pressurized  fractures  than  for  remote  loading,  so 
the  actual  boundary  conditions  imposed  include  zero  displacements  on  the  upper  and 
lower  surfaces,  internal  pressure  on  the  fracture  surface  (p=Eme),  and  the  appropriate 
symmetric  and  interfacial  boundary  conditions.  The  superposed  solutions  are  equivalent 
for  the  constant  shear  stress  interface  model.  Stress  intensity  factors,  remote  loads,  and 
fiber  stresses  have  been  computed  for  three  distinct  crack  configurations  (See  Figure  4.4) 
subject  to  a  range  of  interfacial  strengths,  including  2  MPa,  20  MPa,  and  40  MPa.  This 
range  was  chosen  to  include  both  brittle  and  toughened  failure  modes  and  to  coincide 
with  experimental  data.  Results  of  these  analyses  are  presented  below  in  Section  4.2. 
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Figure  4.3  Matrix  Crack  Initiation  Model;  Matrix,  fiber  and  fracture  discretization  is 
shown  for  the  matrix  crack  intitiation  model  constructed  for  analysis  of  matrix  crack 
initiation  subject  to  varying  interfacial  shear  stress  (2-40  MPa).  The  complete  model 
consists  of  roughly  2500  degrees  of  freedom. 
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4.2  Fracture  Parameter  Results 


4.2.1  Matrix  Cracking  Stress 


The  estimated  matrix  cracking  stresses  for  both  models  are  plotted  in  normalized 
form  in  Figure  4.5  and  tabulated  as  effective  toughness  in  Table  4.2. 


Table  4.2  Effective  Toughness  for  Small  Crack  Propagation; 


Keff/KiC 


c/r\T 

2  MPa 

20  MPa 

40  MPa 

0.5 

5.95 

6.90 

7.30 

1.0 

5.80 

6.76 

7.19 

2.0 

5.42 

6.27 

7.09 

Because  the  toughenig  mechanisms  and  notch  insensitivity  of  fiber-reinforced 
ceramics  rely  on  the  presence  of  fiacture-bridging  fibers,  the  isolation  of  the  fibers  from 
crack-tip  stresses  is  critical  to  successful  use  of  these  materials  in  structural  applications. 
The  axial  fiber  stresses  are  estimated  as  a  function  of  the  shape  function  derivatives  and 
the  collocation  point  displacements  and  computed  for  each  of  the  modeled  fibers. 
Typical  stress  contours  are  shown  in  Figure  4.6  for  the  matrix  crack  initiation  model  and 
the  maximum  is  plotted  as  a  function  of  characteristic  crack  length  for  the  interfacial 
conditions  considered  (See  Figure  4.7). 
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Figure  4.5  Matrix  Cracking  Stress  for  Small  Cracks;  Normalized  values  for  the  matrix 
cracking  stress  are  plotted  along  with  'intermediate'  crack  solutions. 


4.3  Discussion  of  Results 

Recent  results  derived  by  Meda  and  Steif  suggest  that  steady-state  matrix  cracking 
may  occur  in  ceramic  composite  materials,  a  posssibUity  that  holds  great  promise  for 
structural  engineers.  The  results  of  this  investigation  support  the  steady-crack  hypothesis 
and  appear  to  be  consistent  with  these  'intermediate'  crack  length  models  (See  Figure  4.5). 
In  particular,  the  effective  toughening  due  to  crack  pinning  and  interfacial  slip  have  been 
shown  to  be  significantly  higher  than  previous  estimates.  Fiuthermore,  the  effective 
toughness  prior  to  fiber  failure  is  dependent  on  the  interfacial  shear  for  small  crack 
growth.  This  is  an  anticipated  result  of  the  crack  pinning  in  a  region  of  interfacial  sliding 
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Figure  4.7  Peak  Axial  Fiber  Stress  for  Proximal  Fibers;  The  peak  axial  fiber  stress  due 
to  small  matrix  crack  extension  is  plotted  for  proximal  fibers. 


and  of  the  relative  scale  of  the  fracture. 

These  results  show  that  failure  of  the  first  or  second  fiber  (depending  on  the  load 
transfer  and  interface  slip)  can  occur  in  the  presence  of  small  matrix  flaws.  Relatively 
low  interfacial  strengths  are  required  to  safely  isolate  the  fibers  from  the  crack-tip 
stresses.  For  the  SiC/LAS  system  modeled,  interface  sliding  stresses  slightly  above  2 
MPa  can  lead  to  failure  at  room  temperature.  This  is  consistent  with  experimental 
evidence  and  may  have  serious  implications  for  composites  reinforced  with  embrittled 
fibers  with  degraded  strengths  [7].  General  ceramic  composite  systems  are  expected  to 
demonstrate  similar  behavior. 
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5.1  Conclusions 

A  matrix  crack  model  has  been  constructed  for  small  crack  growth  based  on 
experimental  observations  of  crack  location  and  matrix  cracking  stress  in  silicon  carbide 
reinforced  lithium  alumino-silicate  (SiC/LAS)  [4,5].  Kim  and  Pagano  have  observed 
matrix  crack  initiation  in  the  region  between  the  reinforcing  fibers  and  on  the  scale  of  the 
fiber  spacing  [4].  More  importantly,  this  initial  fracture  propagation  occurs  prior  to 
observable  non-linearities  in  the  composite  load-deflection  relationship  which  signify 
either  brittle  failure  or  toughened  fracture. 

This  situation  has  been  analyzed  in  conjunction  with  interface  slip  using 
numerical  techniques  and  more  recently  using  an  axisymmetric  distributed  spring  model 
[18,27].  This  improved  axisymmetric  model  suggests  the  possibility  of  stable  fracture 
growth  for  small  crack  sizes  and  weak  interfacial  conditions.  Stable  matrix  crack 
propagation  is  a  possibility  not  previously  considered  feasible  but  with  significant 
implications  for  the  design  of  toughened  ceramic  materials. 

The  results  of  this  investigation  support  the  steady-crack  hypothesis  and  are 
consistent  with  'intermediate'  crack  length  models  (See  Figure  4.5).  In  particular,  the 
effective  toughening  due  to  crack  pinning  and  interfacial  slip  have  been  shown  to  be 
significantly  higher  than  previous  estimates.  Small  fractures  propagating  among 
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reinforcing  fibers  are  shielded  from  the  applied  loads  as  expected.  Though  it  is  difficult 
to  confirm  with  the  data  collected  to  date,  the  three-dimensional  nature  of  composites  on 
this  scale  appears  to  be  important  and  slows  the  decrease  in  matrix  cracking  stress  as  the 
fracture  propagates.  Despite  this  temporary  change,  the  magnitudes  of  the  critical  stress 
are  close  to  those  derived  by  Meda  and  Steif  using  the  axisymmetric  model  [27], 

The  propagation  behavior  of  small  matrix  flaws  is  critical  to  the  toughening 
mechanisms  available  for  large  fractures  and  is  dependent  on  the  interfacial  strength.  As 
shown  in  Figure  4.7,  failure  of  the  first  or  second  fiber  (depending  on  the  load  transfer 
and  interface  slip)  can  occur  at  relatively  low  interfacial  strengths.  For  the  SiC/LAS  sys¬ 
tem  modeled,  interface  sliding  stresses  above  12  MPa  can  lead  to  failure  at  room 
temperature  and  guarantee  failure  of  degraded  fibers  [7].  General  ceramic  composite 
systems  are  expected  to  demonstrate  similar  behavior. 

On  a  different  level,  this  investigation  has  developed  and  demonstrated  an  effi¬ 
cient  computational  technique  for  analysis  of  three-dimensional  fracture  growth  in  com¬ 
posite  media.  The  technique  appears  to  be  particularly  well-suited  for  microstractural 
models  such  as  the  matrix  crack  model  analyzed  in  this  investigation  and  ready  for 
extension  to  thermo-elastic  fatigue. 

5.2  Recommendations  for  Further  Work 

While  this  research  project  has  resolved  some  issues  regarding  matrix  crack 
initiation  in  fiber-reinforced  ceramics,  it  leaves  many  issues  unresolved  -  by  the  nature  of 
it  finite  scope  -  and  poses  as  many  more.  Some  of  these  issues  are  mentioned  briefly  here 
with  the  hope  that  continued  research  may  soon  resolve  them. 

The  results  presented  are  for  a  single  crack  propagating  from  a  small  matrix  flaw. 
While  it  is  possible  that  a  large,  'steady-state*  crack  can  form  by  extension  of  a  single  flaw 
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,  it  is  more  likely  to  be  form  by  coalescence  of  many  small  flaws.  The  interaction 
between  neighboring  flaws  is  an  important  concept  to  the  overall  behavior  of  the  material. 
In  addition,  the  link  between  these  microscopic  material  events  and  the  macrostructural 
behavior  of  the  material  is  critical  to  the  successful  development  of  these  materials. 

Equally  important  are  the  implications  this  damage  development  process  may 
have  on  useful  service  life  and  maintenance.  Toughened  materials  which  can  be  tested 
in-service  during  regular  maintenance  intervals  will  find  greater  use  and  be  more  reliable. 
Proper  tailoring  of  the  interface  offers  some  opportunities  in  this  direction. 

Many  issues  regarding  chemical  interaction  and  environmental  evolution  of  the 
constituent  have  only  recently  been  addressed.  While  it  is  now  possible  to  tailor 
interfaces  to  provide  toughened  material  behavior  and  failure  modes,  the  same  materials 
suffer  structural  failure  at  high  temperatures  and  in  corrosive  environments  [7]. 

Extension  to  harsher  environments  (e.g.  irradiation)  presents  further  complications. 

Despite  these  unresolved  issues,  continued  interest  in  the  development  and 
analysis  of  these  materials  combined  with  efforts  for  the  economical  manufacture  of 
constituent  materials  promises  bright  futures  for  the  use  of  ceramics  for  high-temperature 
structural  applications 
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Appendix  A 
Fundamental  SIBEH  Relations 


Appendix  A  is  included  as  a  supplement  to  Chapters  2  and  3  which  describe  the 
siuface  integral  and  boundary  element  hybrid  (SIBEH)  method.  This  section,  in 
conjunction  with  Appendix  B,  summarizes  the  fundamental  relations  of  the  technique  and 
describes  its  development  in  detail  sufficient  for  rederivation.  The  equations  presented 
include  the  fundamental  solutions,  element  mapping  functions  and  element  shape 
functions  for  both  methods. 

A.l  Fundamental  Solutions  -  Surface  Integral  Method 

The  fundamental  solutions  on  which  the  surface  integral  method  is  based  can  be 
derived  from  point-force  elasticity  solutions  for  three-dimensions  and  are  analogous  to 
the  dislocation  loop  solutions  useful  in  two-dimensional  fracture  mechanics.  As 
described  in  Chapter  2,  a  combination  of  force  dipoles  (or  multipole)  is  used  to  represent 
infinitesimal  tensile  and  shear  fracture  events. 

These  influence  functions  have  been  developed  from  Kelvin's  solution  for  a 
concentrated  point-force  acting  in  an  infinite  isotropic  domain  and  from  Rongved's 
solutions  for  a  point  force  acting  in  one  of  two  bonded  semi-infinite  regions  [43].  These 
relations  express  the  stress  and  displacements  in  media  siurounding  a  point  force  in  terms 
of  Papkovitch-Neuber  potential  functions,  B  and  p. 
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(A.1) 


u  =  B- 


V(rB+)3) 
4(1 -v) 


(A.2) 


with  ^V^B  =  0  andjuV^^  =  0 


where  a  and  u  represent  the  stresses  and  displacements  at  some  point  in  the  material 
described  by  position  vector  r,  p  and  v  are  material  parameters,  and  Bx,By,  Bz,  and  P  are 
potential  functions  which  depend  upon  the  mechanical  material  properties  and  the 
locations  of  the  sampling  point  and  the  force  application  point. 

A.  1 . 1  Kelvin's  Point-Force  Solution 


For  the  simplest  case  of  a  concentrated  point-force,  P,  acting  in  an  infinite 
homogeneous  medium,  the  potential  functions  are: 


B:  = 


4;rp||r|| 


,i3  =  0 


(A.3) 


A.  1.2  Rongved's  Point-Force  Solution 

For  the  more  complicated  case  of  a  concentrated  point-force  acting  in  one  of  two 
bonded  semi-infinite  domains  (Figure  A.l),  several  sets  of  potential  functions  are 
required  to  completely  describe  the  stress  and  displacement  fields  in  both  materials  [43]. 
For  a  concentrated  force  acting  perpendicular  to  the  bimaterial  interface: 


B,=By=0  and  B>B;  =  0 


(A.4a) 
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Material  1 


nATlAIAL  2 


Figure  A.l  Concentrated  Force  Acting  in  Bimaterial  Domain;  This  figure  depicts  the 
reference  coordinate  frame  used  in  Rongved's  elasticity  solution  for  a  concentrated  point- 
force  acting  in  one  of  two  semi-infinite,  bonded  materials. 
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P'=P. 


(l-V) 


(1-v) 
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(A.4c) 


with  gi  = 


(l-v)c 
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Ri  =  [(x  -  a)^  +  (y  -  bf  +  (z  -  c)^], 
and  Rj  =  [(x  -  a)^  +  (y  -  b)^  +  (z  +  c)^  j 


(A.4d) 


where  the  unprimed  and  primed  terms  pertain  to  the  two  distinct  material  regions  z  >  0 
and  z  <  0  respectively  and  the  parameters  (x,y,z)  and  (a,b,c)  are  the  sampling  and  force- 
application  point  coordinates. 

Similar  potential  functions  have  been  derived  for  a  force  acting  parallel  to  the 
interface  and  are  included  here  for  a  the  case  of  a  concentrated  force  acting  in  the 
direction  of  the  x-axis: 

B  =  ^2- j— + £riL  J-| 
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A.  1.3  Displacement  Discontinuity  Solutions 

From  these  potential  representations,  stress  and  displacement  solutions  can  be 
derived  for  the  force  dipoles  which  comprise  the  displacement  discontinuity  model. 
These  dipole  solutions  are  derived  using  a  limiting  process  for  two  parallel,  opposing 
forces,  P,  separated  by  an  infinitesimal  distance,  h.  Taking  the  limit  of  the  stress  and 
displacement  fields  surrounding  this  force  couple  as  h— >0  (with  the  dipole  magnitude, 
Ph,  held  constant)  gives  the  dipole  solutions  in  terms  of  derivatives  of  the  potential 
functions  taken  with  respect  to  the  distinguishing  distance,  h  [35].  Equations  (A.6)  and 


81 


Appendix  A:  SIBEH  Relations 


(A.7)  represent  the  stress  and  displacement  solutions  for  a  dipole  composed  of  forces  in 
the  K  direction  combined  using  a  limiting  process  in  the  X  direction. 


+0-2v;(b'S  +B'«) 

(A.6) 

(A.7) 

where  i,jdc,m  €  {x,y,z}  and  K,A,e  {a,b,  c} 


Equations  (A.6)  and  (A.7)  are  expressed  in  indicial  notation  so  that  a  comma  in  the 
subscript  denotes  differentiation  with  respect  to  the  trailing  variables  and  repeated  indices 
signify  summation.  These  dipole  relations  are  combined  to  form  multipoles  (as  described 
in  Chapter  2  and  depicted  in  Figure  2.1)  and  are  then  converted  to  displacement 
discontinuity  solutions  with  a  combination  of  material  parameters  [35].  For  a  tensile 
multipole  aligned  with  the  z-axis: 
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For  a  shear  dipole  acting  in  the  x-z  directions: 


(A.10) 
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A  similar  derivation  process  may  be  followed  for  fractures  on  or  near  planar 
bimaterial  interfaces  using  Rongved's  solutions  (Section  A.  1.2).  Once  the  bimaterial 
displacement  discontinuity  solutions  have  been  formed,  interfacial  crack  solutions  can  be 
generated  by  taking  the  limit  as  the  infinitesimal  fracture  event  approaches  the  interface 
from  either  material.  (The  two  distinct  solutions  produce  the  same  stress  and 
displacement  fields  at  this  point.)  Although  somewhat  more  complicated,  these  functions 
can  be  treated  with  the  same  integration  and  solution  schemes  outlined  for  fractures  in 
homogeneous  media. 

Because  of  the  complexity  involved  in  these  derivations,  the  symbolic 
manipulator  MACS  YMA  has  been  used  for  accurate  derivation  of  the  fundamental 
solutions.  However,  their  uniqueness  precludes  straightforward  verification.  For  this 
analysis,  verification  involved  derivation  and  comparison  of  two  distinct  influence 
function  sets,  checks  of  symmetric  properties,  and  comparison  with  related,  previously 
derived  fundamental  solutions  [35]. 


The  primary  boundary  element  influence  functions  for  three-dimensional 
elastostatics  are  also  derived  from  Kelvin's  elasticity  solution  for  a  concentrated  force 
acting  in  an  infinite  isotropic  domain  [51].  Since  complete  derivations  are  presented  in 
many  reference  texts,  only  the  final  function  forms  are  included  here.  For  the  stress  and 
displacement  fields  surrounding  a  unit  point-force  in  the  j-direction: 
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+a-2vj(nirj-njri); 


(A.  13) 


where  u*ij  represents  the  displacement  components  uj  at  point  r  on  the  boundary  element 
surface  induced  by  a  concentrated  unit  force,  Pj,  applied  at  the  origin.  p*ij  represents  the 
corresponding  tractions  at  this  surface  with  a  local  outward  normal  n. 

For  the  boundary  collocation  formulation,  Equations  (A.  12)  and  (A.  13)  are  used 
to  relate  the  boundary  tractions  and  displacements.  As  outlined  in  Chapters  2  and  3,  this 
common  solution  scheme  results  in  a  linear  equation  system  which  can  be  solved  to 
obtain  the  boundary  parameters.  Once  these  boundary  values  have  been  calculated,  the 
derivative  influence  functions  d*  and  s*  can  be  used  to  determine  the  stresses  in  the 
domain  interior  [51].  These  derivative  functions  give  the  stresses  acting  at  an  internal 
point  (e.g.  crack  surface)  as  a  function  of  the  boundary  displacement  and  traction 
distributions  respectively  and  are  derived  by  differentiation  of  the  fundamental  solutions 


u*  and  p*. 
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Note  that  for  these  derivative  relations  the  radial  vector,  r,  extends  from  the  internal  point 
towards  the  external  boundary  and  that  all  derivatives  of  this  vector  are  taken  at  the 
boundary  point 

A.3  Element  Mapping  Functions 


The  element  mapping  functions  incorporated  in  the  matrix  crack  analysis  code  are 
•  based  on  the  bi-quadratic  Lagrange  functions  and  relate  the  local  element  reference  frame 

to  the  global  coordinate  frame  [45]: 


• 

=  o=  1,2,3,. ..,9 

(A.  16) 

where  = 

(A.  17a) 

• 

(A.  17b) 

(A.  17c) 

• 

(A.17d) 

(A.17e) 

• 

(A.17f) 

(A.17g) 

• 

(A.17h) 

(A.17i) 
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Simpler  mapping  functions  are  included  in  this  set  by  condensation.  For  example,  linear 
boundaries  are  approximated  by  parabolically  curved  boundary  segments  for  which  the 
midpoint  is  located  directly  between  the  neighboring  comer  nodal  points.  Similarly, 
triangular  elements  are  included  by  condensing  the  third  and  foiuth  nodal  points  and  by 
combining  the  related  mapping  functions. 

A.4  Element  Shape  Functions 

A  variety  of  local  distribution  functions  have  been  included  for  approximating 
tractions  and  displacements  along  fracture  surfaces  and  component  boundaries.  Three 
shape  functions  have  been  implemented  for  the  surface  integral  model  -  constant, 
constant-linear,  and  crack-tip  elements  [26].  The  corresponding  shape  functions  are 
outlined  below: 


(A.18) 

(A.  19) 

8a)  =  ^,(p(^))8^  +  ^,(pm8'^ 

(A:20) 

with  = 

(A.21) 

where  5^  and  6^  are  the  crack  opening  displacements  at  the  collocation  points  and  Ni(^) 
and  N2(0  are  orthogonal  functions  of  the  actual  crack-tip  radius,  p.  These  functions  are 
based  on  the  first  two  terms  of  the  Williams  expansion  for  fractures  in  homogeneous 
media  and  are  presented  in  complete  form  in  Equation  2. 12. 

For  the  isoparametric  boundary  element  formulation,  the  shape  functions  will  be 
equivalent  to  the  mapping  functions  given  in  section  A.3  above. 


=  for  a=1..9 
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(A.22) 


The  Lagrange  family  of  functions  is  preferred  in  this  case  because  it  results  in  more 
uniform  collocation  point  distributions  and  therefore  greater  solution  accuracy.  Lower 
order  variations  and  other  function  families  (e.g.  serendipity)  can  be  simulated  as  linear 
combinations  of  the  nine-noded  Lagrange  functions  presented  above. 
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Singular  Integration  Schemes 


Appendix  B  has  been  included  as  a  supplement  to  Chapters  2  and  3  to  provide  a 
more  detailed  development  of  the  singular  integration  schemes  used  in  the  SIBEH 
method.  Singular  integral  regularization  procedures  are  outlined  for  both  the  surface 
integral  and  boundary  element  methods.  These  approaches  are  useful  both  for  evaluation 
of  genuinely  singular  integrals  and  for  reduction  of  computational  errors  in  nearly- 
singular  integrals. 


B.l  vSingular  Integration  -  Surface  Integral  Method 


Singular  integration  of  the  fundamental  stress  and  displacement  solutions  for  the 
surface  integral  model  use  similar  approaches  (i.e.  regularization  by  subtraction  of  a  rigid 
body  motion)  but  are  handled  in  slightly  different  manners  due  to  a  difference  in  the 
dominant  singularity  order. 

As  presented  in  Section  2.2,  the  singular  integral  for  stresses  can  be  reduced  to  a 
tractable  form  by  subtracting  a  rigid  body  motion  (i.e.  constant  displacement  over  the 
entire  fracture  plane,  St)  which  contributes  nothing  to  the  stress  fields. 

=  (B.l) 


C„  +  (B.2) 
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using  N;”n£  r‘‘iA=o 


(B.4) 


where  the  original  integral  of  f  (the  fundamental  stress  solution)  is  evaluated  over 
(the  elemental  regions  surrounding  collocation  point  J),  combines  the  element  shape 
functions  associated  with  this  point  J,  is  the  shape  function  value  at  J,  and  n  denotes 


the  normal  direction  associated  with  the  sampling  point. 

The  first  integral  in  Equation  (B.3)  is  now  defined  in  a  Cauchy  principal  value 
sense  and  can  be  evaluated  numerically.  When  the  quadrature  point  coordinates  are  taken 
to  be  symmetric  about  the  singular  point,  Gauss-Legendre  quadrature  orders  as  low  as  2 
can  be  used  to  obtain  accurate  integral  solutions. 

The  remaining  term  in  Equation  (B.3)  involves  integration  of  the  fundamental 
stress  solution  over  the  region  St  (i.c.  the  entire  fracture  plane  excepting  the 
elemental  regions  surrounding  point  J).  Though  impractical  to  evaluate  with  two- 
dimensional  quadrature,  this  integral  is  defined  because  of  the  1/R^  influence  function 
singularity  and  because  of  the  exception  of  point  J  from  the  integration  region.  Further 
simplification  is  attained  by  conversion  to  a  local  polar  coordinate  reference  frame,  which 
reduces  the  dominant  singularity  by  one  order  and  facilitates  analytical  integration  of  the 
radial  terms. 


i.sr<  r’dA = ae 

(B.5) 

• 

=  'Elr,<  /«(■  r.z)!*  as 

(B.6) 

=-'^fr,<e)r^(e,z)de 

(B.7) 

• 

where  rR(0,z)=  [“  y^(r,z)rdr 

•»r(  a) 

(B.8) 
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In  this  case,  the  radial  integral  expressed  in  Equation  (B.8)  is  evaluated  from  the  element 
boundary,  r(0),  to  infinity,  where  the  1/R  singularity  of  r®R(r,z)  reduces  the  contribution 
at  this  bound  to  zero.  To  simplify  the  one-dimensional  integration  in  Equation  (B.7)  for 
parabolically  curved  boundaries,  each  element  boundary  segment  is  represented 


parametrically. 


Evaluation  of  the  fundamental  displacement  solution  for  proximal  points  follows 
a  similar  approach,  though  the  evaluation  of  the  regularizing  integral  is  handled 
differently  because  of  the  lower  order  of  the  dominant  singularity  (l/R^).  Subtraction  of 
a  rigid  body  motion  over  the  entire  fracture  plane  is  not  required.  As  in  the  stress 
solution  case,  a  constant  regularizing  term,  No(J>,  is  removed  to  define  the  Cauchy 
principal  value  integral. 


(B.IO) 


In  this  case,  though,  the  regularizing  integral  is  evaluated  in  a  semi-analytic  fashion  over 
the  elemental  region  surrounding  collocation  point  J. 


r'dA = XJ  de 

fB(d)Ti(e, z)dd 

where  rR('0,z)  =  /iwJ^^  ^  y^(r,z)tds 


(B.ll) 


(B.12) 


(B.13) 
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Existence  of  this  integral  term  requires  that  the  angular  integral  for  a  circular  contour 
(constant  radius)  be  finite. 

limri(e,z^J^j  Yt(B)de]  =  0  (B.14) 

This  property  is  satisfied  by  the  fimdamental  displacement  solutions  for  cases  of  interest 
in  the  surface  integral  and  boundary  element  hybrid  approach. 

B.2  Singular  Integration  -  Boundary  Element  Method 

Two  cases  of  singular  integration  occur  regularly  in  boundary  element  analysis. 
These  involve  evaluation  of  the  primary  fundamental  solutions,  u*  and  p*,  over  an 
elemental  region  coincident  with  the  collocation  point  at  which  tractions  and 
displacements  are  evaluated.  Other  singular  situations  are  possible  in  the  hybrid 
formulation  (e.g.  a  surface  integral  collocation  point  coincident  with  the  boundary 
element  surface)  but  occur  only  in  improperly  posed  problems.  This  section  outlines  an 
integration  scheme  developed  recently  for  general  boundary  element  applications  [44,57]. 

The  implemented  technique  eliminates  the  singularity  of  the  integral  by 
subtracting  a  Taylor  series  expansion  of  the  integrand  about  the  singular  point  For  little 
additional  cost  compared  with  more  traditional  integration  procedures,  the  technique  also 
provides  more  accurate  results  for  proximal  sampling  points.  Although  the  two  methods 
were  developed  independently,  the  singular  integration  method  presented  above  for 
fracture  solutions  can  be  viewed  as  a  specialized  form  of  this  approach.  A  formal 
derivation  and  demonstration  has  been  presented  but  is  outlined  here  for  completeness 
[44].  Because  the  fundamental  solution  singularities  are  one  order  lower  (1/R2  and  1/R 
respectively  for  stresses  and  displacements)  than  those  for  fracture  events,  the  restriction 
to  planar  elements  and  internal  collocation  points  can  be  relaxed. 
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*  Figure  B.l  Local  Coordinate  Frame  used  for  Singular  Integration  Scheme  Local  co¬ 
ordinates,  are  defined  such  that  the  Ci"C2  plane  is  tangent  to  the  element  smface  at  the 
nearest  collocation  point  and  the  origin  is  the  projection  of  the  sanipling  point  into  this 
plane.  The  regularizing  integral  terms  are  evaluated  and  integrated  in  this  tangent  plane. 
(Adapted  from  Reference  [44]) 

• 

To  simplify  the  regularization  procedure  for  the  general  case  of  curved  boundary 
elements,  the  integral  is  evaluated  in  a  coordinate  frame  defined  by  a  plane  (Ci"C2  plane) 
®  tangent  to  the  element  at  the  collocation  point,  Q,  nearest  to  the  source  point,  P.  Further, 

the  coordinate  frame  is  defined  such  that  the  -direction  and  element-defined 
direction  are  parallel  and  the  origin  is  located  at  the  projection  of  the  source  point  onto 
•  this  plane  (See  Figure  B.l).  This  conversion  facilitates  the  development  and  integration 

of  the  regularizing  terms  in  the  local  tangent  plane. 
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Projections  of  points  on  the  element  surface,  or  'image'  points,  can  be  expressed  as 
a  linear  combination  of  mapping  functions  and  the  nodal  point  projections  and  are 
designated  here  with  a  prime: 

0=1 . 9, 1 =1.2  (B.16) 


The  regularizing  terms  used  to  reduce  the  integral  singularity  consist  of  the 
fundamental  solution,  a  first-order  Taylor  series  expansion  of  the  displacements  or 
tractions  (Lq),  and  the  differential  area  dA'  all  evaluated  at  the  image  points  in  the  tangent 
plane. 

4„/Vfp.f>NrfwA = - 

/^.../Vfp.fmffldA-  (B.17) 


where  f*(Co’^)  f*  (Co )  represent  either  of  the  fundamental  solutions  evaluated 

on  the  element  siuface  and  on  the  tangent  plane  (respectively)  and. 


oQi 


(Ct-G) 


(B.18) 


As  a  result  of  the  coordinate  transformation  described  earlier,  the  derivative  terms  in 
Equation  (B.18)  can  be  estimated  at  the  point  Q  as  a  linear  function  of  the  mapping 
function  derivatives  and  nodal  point  coordinates. 

This  expansion  about  Q  can  be  shown  to  completely  negate  the  singularity  in  the 
boundary  element  fundamental  solutions,  resulting  in  a  non-singular  boundary  element 
formulation.  In  addition,  evaluation  of  the  regularizing  terms  within  the  tangent  plane 
allows  them  to  be  expressed  in  an  analytically  integrable  form.  Using  an  approach 
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similar  to  that  developed  above  for  fracture  solutions,  the  planar  function  forms  are 
converted  to  a  local  cylindrical  coordinate  frame  and  separated  into  products  of  radial  and 
angular  terms.  The  radial  terms  are  analytically  integrated  from  the  origin  to  the  element 
boundary  image.  Finally,  the  combined  angular  terms  can  be  integrated  using  low-order 
one-dimensional  quadrature. 

Although  this  integration  scheme  permits  evaluation  of  stress  and  displacement 
fields  nearer  to  the  element  surface  than  previously  possible,  it  is  not  generally  the  most 
cost-effective  approach.  For  source  points  further  from  the  element  than  the  element 
dimension,  straightforward  two-dimensional  numerical  integration  gives  equivalent 
results  at  less  expensive. 


